MASARYKOVA UNIVERZITA Přírodovědecká fakulta Ústav matematiky a statistiky Diplomová práce Brno 2019 Natália Gažová MASARYKOVA UNIVERZITA Přírodovědecká fakulta Ústav matematiky a statistiky Modelování extrémních událostí Diplomová práce Natália Gažová Vedoucí práce: RNDr. Radim Navrátil, PhD. Brno 2019 Bibliografický záznam Autor: Bc. Natália Gažová Přírodovědecká fakulta, Masarykova univerzita Ústav matematiky a statistiky Název práce: Modelování extrémních událostí Studijní program: Matematika Studijní obor: Finanční matematika Vedoucí práce: RNDr. Radim Navrátil, PhD. Akademický rok: 2018/2019 Počet stran: ix + 60 Klíčová slova: teorie extrémních hodnot; metoda blokových maxim; zobecněné rozdělení extrémních hodnot; sféra přitažlivosti; peaks-over-threshold; zobecněné Paretovo rozdelení Bibliografický záznam Autor: Bc. Natália Gažová Prírodovedecká fakulta, Masarykova univerzita Ústav matematiky a štatistiky Názov práce: Modelovanie extrémnych udalostí Študijný program: Matematika Študijný odbor: Finančná matematika Vedúci práce: RNDr. Radim Navrátil, PhD. Akademický rok: 2018/2019 Počet strán: ix + 60 Kľúčové slová: teória extrémnych hodnôt; metóda blokových maxím; zovšeobecnené rozdelenie extrémnych hodnôt; sféra príťažlivosti; peaks-over-threshold; zovšeobecnené Paretovo roz- delenie Bibliographic Entry Author: Bc. Natália Gažová Faculty of Science, Masaryk University Department of Mathematics and Statistics Title of Thesis: Modelling extremal events Degree Programme: Mathematics Field of Study: Financial Mathematics Supervisor: RNDr. Radim Navrátil, PhD. Academic Year: 2018/2019 Number of Pages: ix + 60 Keywords: extreme value theory; block maxima method; generalised extreme value distribution; domain of attraction; peaks-over-threshold; generalised Pareto distribution Abstrakt V této diplomové práci se věnujeme analýze extrémních událostí. V modelování využíváme dvě metody, konkrétně metodu blokových maxim a metodu založenou na překročení určité meze. Úvodní část obsahuje teoretické seznámení s oběma metodami, což je důležité pro dostatečné pochopení při dalších aplikacích. Praktická část se soustředí na analýzu pojistných dat oběma způsoby s následným porovnáním výsledků. Abstrakt V tejto diplomovej práci sa venujeme analýze extrémnych udalostí. V modelovaní využívame dve metódy, konkrétne metódu blokových maxím a metódu založenú na prekročení určitej medze. Úvodná časť práce obsahuje teoretické zoznámenie sa s oboma metódami, čo je dôležité pre dostatočné pochopenie pri ďalších aplikáciach. Praktická časť sa sústredí na analýzu poistných dát oboma spôsobmi s následným porovnaním výsledkov. Abstract This thesis deals with extreme value analysis. Two approaches to modelling extreme events are concerned in this text – block maxima method and peaks-overthreshold. In the first part, a theoretical background is given which is essential for further applications. Afterwards, modelling of insurance data is studied. Both above-mentioned methods as well as the comparison of results are included in the analysis. Poďakovanie Na tomto mieste by som chcela poďakovať vedúcemu tejto diplomovej práce, RNDr. Radimovi Navrátilovi, PhD., za ochotu a svoje cenné pripomienky, ktoré mi počas vedenia práce poskytol a za čas, ktorý mi venoval. Prohlášení Prohlašuji, že jsem svoji diplomovou práci vypracovala samostatně s využitím informačních zdrojů, které jsou v práci citovány. Brno 20. prosince 2018 . . . . . . . . . . . . . . . . . . . . . . . . . . . Natália Gažová Obsah Úvod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Kapitola 1. Teória extrémnych hodnôt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Kapitola 2. Metóda blokových maxím . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Sféra príťažlivosti rozdelení extrémnych hodnôt . . . . . . . . . . . . . . . . 9 2.2 Max-stabilita rozdelení extrémnych hodnôt . . . . . . . . . . . . . . . . . . . 16 2.3 Modelovanie metódou blokových maxím . . . . . . . . . . . . . . . . . . . . . . 18 2.3.1 Odhad parametrov GEV rozdelenia . . . . . . . . . . . . . . . . . . . . . 18 2.3.2 Spracovanie dát . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Kapitola 3. Peaks-over-threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1 Stabilita rozdelení pre modelovanie excesov . . . . . . . . . . . . . . . . . . . 29 3.2 Metóda POT pri praktickom spracovaní dát . . . . . . . . . . . . . . . . . . . 32 3.2.1 Odhad parametrov GP rozdelenia . . . . . . . . . . . . . . . . . . . . . . 32 3.2.2 Spracovanie dát . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Kapitola 4. Analýza dát . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1 Dáta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2 Modelovanie metódou blokových maxím . . . . . . . . . . . . . . . . . . . . . . 39 4.3 Modelovanie POT metódou . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.4 Porovnanie metódy blokových maxím a POT metódy . . . . . . . . . . . . 50 4.5 Parametrický model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Závěr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Seznam použité literatury . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 – viii – Úvod Extrémne udalosti, ako už ich názov napovedá, nie sú spojené s bežnými situáciami, no ich výskyt znamená veľké pohyby, vysoké straty či vysoký dopad v určitej sfére záujmu. Preto je ich štúdium zaujímavé z rôznych hľadísk a pre rôzne skupiny po- pulácie. Typickým príkladom extrémnej udalosti je napr. storočná povodeň, t.j. taká povodeň, ktorá sa vyskytne v dlhodobom priemere jedenkrát za sto rokov. Ďalšou oblasťou bohatou na tieto javy je meteorológia, ktorá zahŕňa pôsobenie hurikánov, rekordných zmien teploty, či vlny horúčav. Práve v oblasti meteorologických a klimatických javov sa analýza extrémnych udalostí začala rozvíjať. Extrémne javy však figurujú aj vo finančníctve, a to napríklad v podobe výraznej zmeny hodnoty aktíva za určitý čas. Veľkou sférou, kde sa metódy taktiež využívajú je poisťovníctvo, a to najmä jeho neživotná časť. V niektorých portfóliách sa môže príležitostne vyskytnúť vysoký poistný nárok, ktorý môže ohroziť jeho stabilitu. Okrem takých udalostí ako hurikány, zemetrasenia či letecké nehody sa vysoké poistné nároky vyskytujú aj v iných oblastiach. Patria medzi ne poistné nároky z automobilového poistenia či z poistenia proti požiarom (priemyselné požiare a podobne). Je preto prínosné sa proti všetkým takýmto okolnostiam brániť vopred, s čím pomáha práve analýza extrémnych udalostí. Prvá kapitola práce poskytuje veľmi stručné zoznámenie s metódami, ktoré sa na modelovanie extrémnych situácii využívajú. Prvá z nich – metóda blokových maxím je podrobnejšie študovaná v kapitole 2. Okrem iného popisuje aj sféry príťažlivosti a overuje max-stabilitu jednotlivých rozdelení extrémnych hodnôt. Ďalej je v nej pozornosť venovaná odhadu parametrov rozdelenia extrémnych hodnôt a riešeniu praktických úloh spojených s analýzou dát, ako napr. hľadanie doby návratu extrémneho javu. Nasledujúca kapitola 3 potom popisuje podobné charakteristiky pre metódu peaks-over-threshold. V poslednej kapitole 4 je priestor venovaný analýze dát z oblasti poisťovníctva. Aplikované sú obe zmienené metódy, ktoré sú na záver porovnané. – ix – Kapitola 1 Teória extrémnych hodnôt Klasická štatistika sa, zjednodušene povedané, zameriava na „priemerné správanie nejakého náhodného procesu, teda na javy realizujúce sa v okolí strednej hodnoty nejakého rozdelenia. Na druhú stranu, teória extrémnych hodnôt (ďalej len EVT, z angl. Extreme Value Theory) sa sústredí na extrémne, ojedinelé a vzácne javy, ktoré sa spravidla realizujú na chvostoch rozdelení. Z pohľadu štatistiky ide o tzv. odľahlé pozorovania. Ukazuje sa totiž, že pravdepodobnosť nastania extrémnej udalosti je v skutočnosti vyššia, ako sa predpokladá na základe teoretického rozdelenia pravdepodobnosti, t.j. že rozdelenia majú ťažšie chvosty. Cieľom EVT je skúmanie správania týchto odľahlých pozorovaní. Zaujímať nás napríklad môže stanovenie veľkosti extrémneho javu, ktorý nastane raz za zvolenú frekvenciu, prípadne naopak, stanovenie priemernej frekvencie výskytu extrémneho javu danej veľkosti a ďalšie. To znamená, že musíme nájsť vhodnú matematickú metódu, aby sme mohli vysvetliť okolnosti, ktoré sa vyskytujú s relatívne malou pravdepodobnosťou, no majú vysoký dopad. Najčastejšie sa na modelovanie týchto situácií využívajú dve metódy, a to metóda blokových maxím a metóda peaks-over-threshold (POT). Metóda blokových maxím Podstatou metódy blokových maxím je rozdelenie získaných pozorovaní do určitého počtu blokov, pričom v každom bloku zistíme maximálne pozorovanie. Ideou potom je nemodelovať samotné pozorovania, ale získané blokové maximá. Obrázok 1.1 ilustruje koncept metódy blokových maxím. Pozorovania sú v tomto prípade rozdelené do štyroch blokov po troch pozorovaniach. Prístup k modelovaniu dát pomocou blokových maxím má však jednu zjavnú nevýhodu - v rámci bloku uvažujeme len jednu extrémnu hodnotu. Vo všeobecnosti môže nastať situácia, kedy druhá najvyššia hodnota daného bloku je vyššia ako maximum nejakého iného, čím v rámci analýzy strácame cenné extrémne dáta. – 1 – Kapitola 1. Teória extrémnych hodnôt 2 q q q q 1 2 3 4 Obr. 1.1: Blokové maximá. Obr. 1.2: Peaks-over-threshold Peaks-over-threshold Spomenutý nedostatok metódy blokových maxím odstraňuje metóda peaks-over-threshold (POT). Je považovaná za modernejší spôsob analýzy extrémnych udalostí a je založená na prekročení určitej pevne zvolenej medze. Opäť sa nebudeme snažiť modelovať dané pozorovania. Namiesto toho sa budeme zaoberať prekročeniami pozorovaní (excesmi) vopred danej hranice, v prípade, že k tomuto prekročeniu došlo. Tým sú dáta využité efektívnejšie ako v predchádzajúcom prípade. Na obrázku 1.2 je táto hranica vyobrazená prerušovanou čiarou. V nasledujúcich kapitolách si obe metódy rozoberieme podrobnejšie. Zistené poznatky ďalej využijeme pri spracovaní dát. Kapitola 2 Metóda blokových maxím Majme postupnosť náhodných veličín X1, X2, . . . , Xn, pričom jednotlivé náhodné veličiny sú nezávislé a rovnako rozdelené s distribučnou funkciou F (ďalej ich budeme skrátene označovať ako i.i.d. náhodné veličiny, z angl. independent identically distributed). Označme maximum tohto náhodného výberu Mn = max(X1, . . . , Xn), pričom toto maximum chceme ďalej skúmať. Obdobne by sme mohli vyšetrovať vlastnosti minima, pričom je v takom prípade potrebné previesť transformáciu min(X1, . . . , Xn) = − max(−X1, . . . , −Xn). Zaujíma nás, ako toto maximum modelovať, aké rozdelenie je vhodné využiť. Preto s využitím i.i.d. vlastnosti náhodných veličín X1, . . . , Xn zostrojme distribučnú funkciu pre Mn: Gn(x) = P(Mn ≤ x) = P(X1 ≤ x, . . . , Xn ≤ x) = n i=1 P(Xi ≤ x) = Fn (x). (2.1) Intuitívne očakávame, že maximálne extrémne situácie sa odohrávajú na pravom konci distribučnej funkcie F, čo nám potvrdzuje aj nasledujúca veta. Veta 2.0.1. Nech Mn je maximum i.i.d. náhodných veličín X1, . . . , Xn a xF je pravý koncový bod rozdelenia s distribučnou funkciou F definovaný ako xF := sup{x ∈ R, F(x) < 1}. Potom pre n → ∞ platí Mn s.i. −→ xF . Dôkaz. Musíme rozlíšiť dve situácie. Pre x < xF , x ∈ R máme lim n→∞ P(Mn ≤ x) = lim n→∞ Fn (x) = 0, – 3 – Kapitola 2. Metóda blokových maxím 4 a zároveň pre x ≥ xF a xF < ∞ platí P(Mn ≤ x) = Fn (x) = 1. Z toho okamžite plynie konvergencia maxima Mn k pravému koncovému bodu xF v pravdepodobnosti. Nakoľko však postupnosť maxím Mn pre n → ∞ je neklesajúca, môžeme povedať, že Mn konverguje k xF skoro iste. Z predchádzajúcej vety môžeme vidieť, že limitné rozdelenie maxima v sebe nenesie žiadnu informáciu, nakoľko má degenerovanú distribučnú funkciu. Ponúka sa možnosť normovať maximá vhodnými postupnosťami konštánt. Princíp normovania je podobný tomu, ktorý využíva centrálna limitná veta. Hľadajme preto nedegenerovanú distribučnú funkciu, označme ju H(x), ktorá bude limitnou distribučnou funkciou normovanej náhodnej veličiny Mn, čo znamená, že nás zaujíma taká H(x), pre ktorú pre každé x ∈ R platí P Mn − bn an ≤ x n→∞ −→ H(x), kde an > 0, bn ∈ R sú vhodne zvolené konštanty. Potom distribučná funkcia normovaných maxím je tvaru P Mn − bn an ≤ x = P(Mn ≤ anx + bn) = Gn(anx + bn) = [F(anx + bn)]n , (2.2) kde Gn(anx + bn) je distribučná funkcia nenormovaných maxím v danom bode. Ak konštanty an > 0, bn ∈ R existujú, potom je H(x) distribučná funkcia tzv. rozdelenia extrémnych hodnôt, o čom hovorí nasledujúca veta. Veta 2.0.2 (Fisherova-Tippettova). Nech X1, . . . , Xn sú i.i.d. náhodné veličiny, Mn = max{X1, . . . , Xn}. Ak existujú konštanty an > 0, bn ∈ R a nedegenerovaná distribučná funkcia H(x) taká, že P Mn − bn an ≤ x n→∞ −→ H(x), potom má H(x) jeden z nasledujúcich tvarov: G0(x) = exp{−e−x }, x ∈ R, G1(x) = 0, x ≤ 0 exp{−x−α }, x > 0, G2(x) = exp{−(−x)α }, x ≤ 0, 1, x > 0, pričom α > 0 sa nazýva parameter tvaru. Kapitola 2. Metóda blokových maxím 5 −5 −4 −3 −2 −1 0 0.00.51.01.5 Weibullovo rozdelenie x hustota α 1 3 5 −4 −2 0 2 4 0.00.51.01.5 Gumbelovo rozdelenie x hustota 0 1 2 3 4 5 0.00.51.01.5 Frechetovo rozdelenie x hustota α 1 3 5 Obr. 2.1: Hustoty rozdelení extrémnych hodnôt v štandardizovanom tvare. Dôkaz. Dôkaz vety je pomerne náročný, preto ho na tomto mieste nebudeme uvádzať. Náčrt dôkazu je možné nájsť napr. v Embrechts a kol. [7] v kapitole 3. Celý dôkaz v úplnej všeobecnosti prvýkrát uviedol Gnedenko v roku 1943, pozri [10]. Fisherova–Tippetova veta teda inými slovami vraví, že limitné rozdelenie normovaných maxím nejakého náhodného výberu musí byť jedného z troch zmienených typov, nezávisle na pôvodnom rozdelení náhodného výberu. Rozdelenie s distribučnou funkciou G0(x) sa nazýva Gumbelovo rozdelenie, rozdelenie s distribučnou funkciou G1(x) je Fréchetovo rozdelenie, pričom tieto sa využívajú pri modelovaní maxím. Rozdelenie s distribučnou funkciou G2(x) sa nazýva a Weibullovo1 rozdelenie a slúži na modelovanie miním. Tieto tri rozdelenia nazývame rozdelenia extrémnych hodnôt alebo aj EV rozdelenia (Extreme Value). Na obrázku 2.1 sú zobrazené hustoty rozdelení extrémnych hodnôt pre rôzne hodnoty parametru α. Predchádzajúce definície distribučných funkcií rozdelení extrémnych hodnôt sú v štandardizovanom tvare, čo znamená, že parameter polohy µ = 0 a parameter merítka σ = 1. Úplné triedy týchto rozdelení získame zavedením týchto parametrov polohy a merítka, t.j. Gi,µ,σ = Gi x − µ σ , i = 0, 1, 2, µ ∈ R, σ > 0. Jenkinson [14] ukázal, že Gumbelovo, Fréchetovo a Weibullovo rozdelenie môžeme zapísať aj úspornejšie, pomocou jednej distribučnej funkcie zovšeobecneného rozdelenia extrémnych hodnôt. 1 V niektorej literatúre sa nami uvažované Weibullovo rozdelenie používané na modelovanie extrémnych situácií označuje aj ako inverzné alebo extremálne Weibullovo rozdelenie. Kapitola 2. Metóda blokových maxím 6 −4 −2 0 2 4 0.00.20.40.60.81.0 GEV rozdelenie x hustota γ −1 0 1 Obr. 2.2: Hustota zovšeobecneného rozdelenia extrémnych hodnôt. Definícia 2.0.1. Zovšeobecnené rozdelenie extrémnych hodnôt (GEV - generalised extreme value) s parametrom γ ∈ R má distribučnú funkciu v tvare Gγ(x) = exp −(1 + γx)− 1 γ , 1 + γx > 0, γ = 0, (2.3) pričom γ sa nazýva parameter tvaru. Vzťah (2.3) uvádza tvar distribučnej funkcie GEV rozdelenia v štandardizovanom tvare. Opäť môžeme pridať parametre polohy µ ∈ R a merítka σ > 0, pričom platí Gγ,µ,σ(x) = Gγ x − µ σ , 1 + γ(x − µ) σ > 0. Hustota pravdepodobnosti GEV rozdelenia pre rôzne hodnoty parametra γ je zobrazená na obr. 2.2. Ľahko sa môžeme presvedčiť, že definícia GEV rozdelenia je korektná. Po spojitom dodefinovaní hodnoty parametru tvaru γ → 0, a teda následným limitným prechodom získame Gumbelovo rozdelenie, konkrétne lim γ→0 Gγ(x) = lim γ→0 exp −(1 + γx)− 1 γ = = exp − lim γ→0 (1 + γx)− 1 γ = exp −e−x , čo je naozaj distribučná funkcia Gumbelovho rozdelenia. Ďalej pre γ > 0 a γ = 1/α platí Gγ(x) = G1/α(x) = exp − 1 + x α −α , Kapitola 2. Metóda blokových maxím 7 čo je distribučná funkcia Fréchetovho rozdelenia v bode 1 + x/α, kde α > 0. Nakoniec, pre γ < 0 a γ = −1/α dostávame Gγ(x) = G−1/α(x) = exp − 1 − x α α , čím potvrdzujeme, že sa jedná o distribučnú funkciu Weibullovho rozdelenia v bode −1+x/α. GEV rozdelenie teda zahŕňa všetky prípady rozdelení extrémnych hodnôt. Príklad 2.0.1. Majme náhodný výber X1, . . . , Xn z exponenciálneho rozdelenia s parametrom λ = 1. K akému rozdeleniu konverguje maximum tohoto náhodného výberu? Exponenciálne rozdelenie s parametrom λ = 1 má distribučnú funkciu v tvare F(x) = 1 − e−x , x ≥ 0. Konštanty zvolíme v tvare an = 1, bn = ln n. Prečo sme zvolili práve takýto tvar normujúcich konštánt si vysvetlíme v nasledujúcej podkapitole 2.1. Dostávame P Mn − bn an ≤ x = P (Mn − ln n ≤ x) = P (Mn ≤ x + ln n) . To je distribučná funkcia normovaných maxím v bode x + ln n. Aplikovaním vzťahu (2.2) získavame P (Mn ≤ x + ln n) = [F (x + ln n)]n = 1 − e−(x+ln n) n = 1 − 1 n e−x n . Ďalej platí, že 1 − 1 n e−x n → exp −e−x pre n → ∞, čo je distribučná funkcia Gumbelovho rozdelenia (v štandardizovanom tvare), a teda limitným rozdelením maxím exponenciálneho rozdelenia pre vhodne zvolené normovacie konštanty je práve Gumbelovo rozdelenie. Na obrázku 2.3 je táto konvergencia znázornená pre rôzne rozsahy náhodného výberu pri zvolených konštantách an = 1 a bn = ln n. Vidíme, že konvergencia rastie pomerne rýchlo, pričom už rozdelenie maxima výberu s rozsahom n = 20 je veľmi blízke Gumbelovmu rozdeleniu. Na príklade sme si ukázali platnosť Fisherovej-Tippetovej vety pre konkrétne rozdelenie, a teda, že limitné rozdelenie maxima pre vhodné konštanty an > 0, bn ∈ R konverguje k nejakému rozdeleniu extrémnych hodnôt (v prípade, že takéto konštanty existujú). Poznamenajme, že v prípade zvolenia nejakej inej konštantnej hodnoty pre an v príklade 2.0.1, by sa konvergencia ku Gumbelovmu rozdeleniu nezmenila, zmenila Kapitola 2. Metóda blokových maxím 8 −2 −1 0 1 2 0.00.20.40.60.81.0 x Gumbel n = 1 −2 −1 0 1 2 0.00.20.40.60.81.0 x Gumbel n = 5 −2 −1 0 1 2 0.00.20.40.60.81.0 x Gumbel n = 10 −2 −1 0 1 2 0.00.20.40.60.81.0 x Gumbel n = 20 Obr. 2.3: Distribučná funkcia Gumbelovho rozdelenia a normalizovaných maxím náhodného výberu z exponenciálneho rozdelenia pre rôzny rozsah výberu n. Kapitola 2. Metóda blokových maxím 9 by sa však hodnota parametrov rozdelenia. Z toho vyplýva nejednoznačnosť určenia konštánt an > 0, bn ∈ R, pričom so zmenou výberu konštánt sa môže zmeniť aj rýchlosť konvergencie limitného rozdelenia maxím k nejakému z EV rozdelení, pozri [4]. Môžeme preto nájsť viaceré hodnoty, pre ktoré budú normované maximá konvergovať. Výsledné limitné rozdelenia sa však budú líšiť len v parametroch polohy a merítka. 2.1 Sféra príťažlivosti rozdelení extrémnych hod- nôt Ak vhodne štandardizované maximum i.i.d. náhodných veličín konverguje v distribúcii, potom je otázka, ako spoznať, ku ktorému z troch rozdelení extrémnych hodnôt sa tak deje. Definujme najskôr pojem sféry príťažlivosti. Definícia 2.1.1. Nech X1, . . . , Xn je náhodný výber z rozdelenia s distribučnou funkciou F(x). Pokiaľ distribučná funkcia F(x) spĺňa podmienky Fisherovej-Tippetovej vety, t.j. ak pre vhodné konštanty an > 0 a bn ∈ R a pre všetky x ∈ R platí [F(anx + bn)]n → Gγ(x), n → ∞ potom hovoríme, že F(x) patrí do sféry príťažlivosti rozdelenia extrémnych hodnôt s distribučnou funkciou Gγ(x), píšeme F ∈ MDA(Gγ). Nakoľko sa v prípade teórie extrémnych hodnôt zaujímame o rozdelenie maxima, je prirodzené sa zaoberať pravým chvostom rozdelenia náhodného výberu. Preto by sféra príťažlivosti mala závisieť len na tvare a vlastnostiach chvosta rozdelenia a nie na jeho zvyšku. Toto potvrdzuje aj nasledujúca veta. Pred jej uvedením zaveďme označenie F(x) = 1 − F(x) = P(X > x), x ∈ R, pričom funkcia F(x) sa bežne nazýva funkcia prežívania alebo distribučná funkcia chvostu rozdelenia. Veta 2.1.1. Rozdelenie pravdepodobnosti s distribučnou funkciou F(x) patrí do sféry príťažlivosti EV rozdelenia s distribučnou funkciou Gi(x), i = 0, 1, 2, a normovacími konštantami an > 0 a bn ∈ R práve vtedy, ak platí lim n→∞ nF(anx + bn) = − ln Gi(x), x ∈ R. Dôkaz. Podľa definície 2.1.1 a Fisherovej-Tippetovej vety 2.0.2 platí F ∈ MDA(Gi) ⇔ lim n→∞ [F(anx + bn)]n = Gi(x) = 0. Rovnicu môžeme pomocou funkcie prežívania upraviť do tvaru lim n→∞ 1 − F(anx + bn) n = Gi(x). Kapitola 2. Metóda blokových maxím 10 Následne zlogaritmovaním rovnice a využitím rovnosti − ln(1 − x) ≈ x pre x → 0 (túto aproximáciu môžeme využiť, nakoľko limn→∞ F(anx + bn) = 0) získavame lim n→∞ −n ln 1 − F(anx + bn) ≈ lim n→∞ −nF(anx + bn) = ln Gi(x), z čoho vyplýva vzťah z tvrdenia vo vete. Vzhľadom na to, že predpokladáme, že sféra príťažlivosti úzko súvisí s chvostami rozdelení, definujme, že rozdelenie s distribučnou funkciou FX a rozdelenie s distribučnou funkciou FY sú chvostovo ekvivalentné ak pre príslušné funkcie prežívania platí lim x→∞ FX(x) FY (x) = c, kde 0 < c < ∞ je nejaká konštanta. Potom platí, že pokiaľ sú dve rozdelenia chvostovo ekvivalentné, tak patria do rovnakej sféry príťažlivosti, nakoľko konštantu c je možné zahrnúť do normujúcich konštánt an a bn (pozri [15]). Ďalej si popíšeme sféry príťažlivosti pre jednotlivé EV rozdelenia. Fréchetovo rozdelenie Sféra príťažlivosti Fréchetovho rozdelenia má relatívne elegantné vyjadrenie pomocou pomaly sa meniacej funkcie v nekonečne. Preto tento pojem najskôr definujme. Definícia 2.1.2. Kladná Lebesgueovsky merateľná funkcia L na intervale (0, ∞) je • pomaly sa meniaca v nekonečne, ak platí lim x→∞ L(tx) L(x) = 1, t > 0, • pravidelne sa meniaca v nekonečne s indexom ρ ∈ R, ak platí lim x→∞ L(tx) L(x) = tρ , t > 0. Typickým príkladom pomaly sa meniacej funkcie (vynechaním triviálneho prípadu konštantnej funkcie) je napríklad logaritmus ln(x). Môžeme si všimnúť, že regulárne meniacu sa funkciu možno vyjadriť ako súčin mocninnej a pomaly sa meniacej fun- kcie. Podľa vyššie uvedeného predpokladu, sféra príťažlivosti súvisí s chvostami rozdelení. Chvost distribučnej funkcie Fréchetovho rozdelenia má pre α > 0 s využitím Taylorovho rozvoja exponenciály približne tvar Kapitola 2. Metóda blokových maxím 11 G1(x) = 1 − exp −x−α ≈ x−α , x > 0. Preto môžeme približne očakávať, že rozdelenia patriace do MDA(G1) budú mať tzv. ťažké chvosty, ktoré klesajú mocninne. Ako sa od tejto predstavy môžeme „oddialiť formuluje nasledovné tvrdenie. Veta 2.1.2. Rozdelenie s distribučnou funkciou F(x) patrí do sféry príťažlivosti Fréchetovho rozdelenia s distribučnou funkciou G1(x), práve keď platí F(x) = x−α L(x), α > 0 (2.4) pre nejakú pomaly sa meniacu funkciu L. Naviac, ak F ∈ MDA(G1), potom Fn (anx) → G1(x), n → ∞ a normujúce konštanty môžeme voliť v tvare an = F−1 1 − n−1 , bn = 0. (2.5) Dôkaz. Pozri [7], kap. 3. Veta nám hovorí, že sféru príťažlivosti Fréchetovho rozdelenia tvoria distribučné funkcie, ktorých chvosty sú regulárne sa meniace v nekonečne s indexom −α. Nakoľko sme spomenuli, že do MDA Fréchetovho rozdelenia patria rozdelenia s ťažkými chvostami, z hľadiska finančných aplikácií sa táto trieda využíva napr. pri modelovaní trhov s vysokou volatilitou. Príklad 2.1.1. Ukážte, že Paretovo rozdelenie s distribučnou funkciou pre α > 0 F(x) = 1 − x−α , x ≥ 1 patrí do sféry príťažlivosti Fréchetovho rozdelenia. Túto skutočnosť môžeme ukázať dvomi spôsobmi. Distribučná funkcia chvosta Paretovho rozdelenia je tvaru F(x) = x−α , čo môžeme triviálne vyjadriť ako súčin konštantnej funkcie L(x) ≡ 1 a mocninnej funkcie x−α , preto platí (2.4) a F ∈ MDA(G1). Druhým spôsobom je použiť priamy výpočet konvergencie pri zvolených konštantách podľa (2.5). Vyjadrime preto tvar konštanty an pomocou kvantilovej funkcie F−1 (1 − n−1 ) nasledovne: F(x) = 1 − n−1 1 − x−α = 1 − n−1 x = n1/α Kapitola 2. Metóda blokových maxím 12 Zvoľme preto an = n1/α , bn = 0 a overme, či platí konvergencia podľa vety 2.1.1. Konkrétne [F (anx + bn)]n = F n1/α x n = 1 − n1/α x −α n = 1 − x−α n n n→∞ −→ exp −x−α , čo je distribučná funkcia Fréchetovho rozdelenia, čím sme potvrdili, že F ∈ MDA(G1). Do sféry príťažlivosti Fréchetovho rozdelenia patria mnohé ďalšie rozdelenia, pozri [2], [7] alebo [17]. Sú to napr. zovšeobecnené Paretovo, F rozdelenie, ďalej Studentovo t rozdelenie, inverzné gama alebo Cauchyho či logaritmické gama rozdelenie. Weibullovo rozdelenie Fréchetovo a Weibullovo rozdelenie majú medzi sebou úzky vzťah, nakoľko pre α > 0 platí G2(−x−1 ) = G1(x), x > 0. (2.6) Skutočne, G2(−x−1 ) = exp − −(−x−1 ) α = exp −x−α = G1(x). Toto vzájomné vyjadrenie zohráva veľkú úlohu pri odvodení a popísaní sféry príťažlivosti Weibullovho rozdelenia. Veta 2.1.3. Rozdelenie s distribučnou funkciou F(x) patrí do sféry príťažlivosti Weibullovho rozdelenia s distribučnou funkciou G2(x), práve keď je pravý koncový bod rozdelenia xF < ∞ a zároveň F(xF − x−1 ) = x−α L(x), α > 0 (2.7) pre nejakú pomaly sa meniacu funkciu L. Naviac, ak F ∈ MDA(G2), potom Fn (anx + xF ) → G2(x), n → ∞ a normujúce konštanty môžeme voliť v tvare an = xF − F−1 1 − n−1 , bn = xF . (2.8) Dôkaz. Dôkaz nutnej podmienky je náročný, je možné ho nájsť v [19]. Postačujúcu podmienku dokážeme nasledovne. Predpokladajme, že xF < ∞ a platí (2.7). Označme F∗(x) = F(xF − x−1 ), x > 0. Kapitola 2. Metóda blokových maxím 13 Podľa Vety 2.1.2 potom F∗(x) ∈ MDA(G1), pričom normujúce konštanty môžeme zvoliť ako a∗ n = F−1 ∗ (1 − n−1 ) a bn = 0. To znamená, že keďže F∗ patrí do sféry príťažlivosti Fréchetovho rozdelenia, platí Fn ∗ (a∗ nx) → G1(x), a preto využitím vzťahu (2.6) získavame Fn (xF − (a∗ nx)−1 ) → G2(−x−1 ). Nakoniec, substitúciou x = −y−1 získavame Fn xF + y a∗ n → G2(y), y < 0, (2.9) čo je distribučná funkcia Weibullovho rozdelenia v bode y. Zostáva určiť tvar normujúcich konštánt. Z (2.9) okamžite vidíme, že bn = xF . Ďalej an = (a∗ n)−1 , pričom podľa (2.5) je a∗ n = F−1 ∗ (1 − n−1 ) = inf x ∈ R : F∗ (x) ≥ 1 − n−1 . Zavedením substitúcie u = xF − x−1 a ďalšími úpravami získavame a∗ n = inf (xF − u)−1 : F(u) ≥ 1 − n−1 = = xF − inf u : F(u) ≥ 1 − n−1 −1 = = xF − F−1 (1 − n−1 ) −1 , z čoho vyplýva, že an = (a∗ n)−1 = xF − F−1 (1 − n−1 ), čo sme chceli dokázať. Z prechádzajúcej vety vidíme, že rozdelenia patriace do sféry príťažlivosti Weibullovho rozdelenia majú obmedzený pravý koncový bod. Preto sa toto rozdelenie vo veľkej miere nevyužíva na modelovanie situácií vo finančných aplikáciach. Napriek tomu môžu byť rozdelenia tejto sféry vhodné pri nejakých špecifických stratégiách ako napr. zaistenie (hedging), kedy sa určitou kombináciou aktív alebo derivátov dokáže vytvoriť situácia, v ktorej je strata (prípadne zisk) investora limitovaný. Príklad 2.1.2. Rovnomerné spojité rozdelenie na intervale (0, 1) patrí do MDA(G2). Distribučná funkcia je tvaru F(x) = x, pre x ∈ (0, 1). Zrejme potom xF = 1 < ∞. Ďalej overme, či platí (2.7): F(xF − x−1 ) = F(1 − x−1 ) = 1 − (1 − x−1 ) = x−1 , takže podmienka je splnená pre L(x) ≡ 1 a α = 1, z čoho vyplýva, že F ∈ MDA(G2). Pri druhom spôsobe vypočítame normujúce konštanty podľa (2.8), teda bn = xF = 1. Kapitola 2. Metóda blokových maxím 14 Ďalej F(x) = 1 − n−1 , x = 1 − n−1 , a preto an = xF − F−1 (1 − n−1 ) = 1 − (1 − n−1 ) = n−1 . Overíme konvergenciu Fn (anx + bn) = 1 + x n n n→∞ −→ exp(x). Pritom exp(x) je distribučná funkcia Weibullovho rozdelenia pre hodnotu α = 1, a teda F ∈ MDA(G2). Okrem rovnomerného spojitého rozdelenia patrí do MDA(G2) napríklad aj beta rozdelenie, pozri [7]. Gumbelovo rozdelenie Popísať sféru príťažlivosti Gumbelovho rozdelenia je zložitejšie ako v predchádzajúcich dvoch prípadoch. Ak sa pozrieme na chvost jeho distribučnej funkcie, zistíme, že (využitím aproximácie exponenciály prvými dvoma členmi príslušného Taylorovho radu) G0(x) ≈ exp(−x), x → ∞. Preto môžeme očakávať, že chvosty rozdelení, ktoré budú patriť do MDA(G0), budú klesať exponenciálne. Opäť musíme rozhodnúť, koľko môžeme z tejto predstavy zľaviť, aby sa neporušili podmienky sféry príťažlivosti. Do MDA(G0) patrí široká trieda rozdelení, ktoré nemajú rovnaké charakteristiky chvostov. Patrí sem jednak normálne rozdelenie s ľahkými chvostami, ale aj napr. lognormálne rozdelenie s chvostami relatívne ťažkými. Túto triedu popisujeme pomocou tzv. von Misesovych funkcií. Definícia 2.1.3. Nech F je distribučná funkcia náhodnej veličiny X s pravým koncovým bodom xF ≤ ∞. Predpokladajme, že existuje z ∈ R také, že pre F platí F(x) = K exp − x z 1 d(t) dt , z < x < xF (2.10) pre nejakú konštantu K > 0 a funkciu d(x), ktorá je kladná a absolútne spojitá s deriváciou d (x), pre ktorú platí lim x↑xF d (x) = 0. Potom funkcia F sa nazýva von Misesova funkcia. Kapitola 2. Metóda blokových maxím 15 Príkladom von Misesovej funkcie je napr. distribučná funkcia exponenciálneho rozdelenia F s funkciou d(x) = λ−1 , nakoľko F(x) = exp(−λx), x ≥ 0, a tým pádom platí (2.10). Postačujúcou podmienkou náležitosti rozdelenia do Gumbelovej sféry príťažlivosti je, že ak distribučná funkcia F je von Misesova funkcia, potom F ∈ MDA(G0) (pozri [7]). To však nie je jej úplná charakterizácia, kedy je potrebné von Misesovu funkciu určitým spôsobom zovšeobecniť. Veta 2.1.4. Rozdelenie s distribučnou funkciou F patrí do sféry príťažlivosti Gumbelovho rozdelenia práve vtedy, keď existuje nejaké z < xF tak, že F(x) môžeme vyjadriť ako F(x) = K(x) exp − x z g(t) d(t) dt , z < x < xF , (2.11) kde K a g sú merateľné funkcie spĺňajúce lim x↑xF K(x) = K > 0 a lim x↑xF g(x) = 1 a d(x) je kladná a absolútne spojitá funkcia, pre ktorú platí lim x↑xF d (x) = 0. Naviac, ak F ∈ MDA(G0), tak pre n → ∞ Fn (anx + bn) → G0(x) a normujúce konštanty môžeme voliť an = d(bn), bn = F 1 − n−1 , (2.12) pričom možný tvar funkcie d(x) je d(x) = xF x F(t) F(x) dt. Dôkaz. Pozri [19]. Príklad 2.1.3. Ako sme už ukázali v Príklade 2.0.1, exponenciálne rozdelenie s parametrom λ = 1 patrí do sféry príťažlivosti Gumbelovho rozdelenia. Táto skutočnosť sa vzťahuje na exponenciálne rozdelenie pre všetky hodnoty parametru λ > 0. Zvoľme K(x) = 1, g(x) = 1, z = 0 a d(x) = λ−1 . Následne dosadením do vzťahu (2.11) získavame F(x) = exp − x 0 λ dt = exp {−λx} , čo je distribučná funkcia exponenciálneho rozdelenia, a teda F spĺňa danú podmienku a patrí do sféry príťažlivosti Gumbelovho rozdelenia. Kapitola 2. Metóda blokových maxím 16 Náležitosť do MDA(G0) môžeme overiť aj priamym výpočtom pomocou normujúcich konštánt vyjadrených na základe vzťahu (2.12). Takže F(x) = 1 − n−1 , 1 − e−λx = 1 − n−1 , x = ln n λ , preto zvolíme bn = λ−1 ln n a an = d ln n λ = 1 λ . Ďalej stačí overiť, že Fn (anx + bn) konverguje pre n → ∞ ku G0(x), takže Fn (anx + bn) = [1 − exp {−(x + ln n)}]n = 1 − e−x n n pričom 1 − e−x n n n→∞ −→ exp e−x , a preto F ∈ MDA(G0). Do sféry príťažlivosti Gumbelovho rozdelenia patrí okrem exponenciálneho široké spektrum často využívaných pravdepodobnostných rozdelení - normálne, gama, (klasické) Weibullovo, lognormálne, či χ2 rozdelenie. Podobnejšie je táto skutočnosť rozoberaná napr. v [2], [4] alebo [7]. 2.2 Max-stabilita rozdelení extrémnych hodnôt Rozdelenia extrémnych hodnôt majú zaujímavú vlastnosť – sú tzv. max-stabilné. Definícia 2.2.1. Rozdelenie náhodného výberu X1, . . . , Xn sa nazýva max-stabilné, ak pre jeho nedegenerovanú distribučnú funkciu F(x) platí, že existuje vhodná postupnosť konštánt an > 0 a bn ∈ R takých, že [F(anx + bn)]n = F(x) pre všetky x ∈ R. Pripomeňme, že [F(anx + bn)]n z predchádzajúcej definície je distribučná funkcia normovaného maxima náhodného výberu. Vlastnosť max-stability majú všetky rozdelenia extrémnych hodnôt, čo si postupne ukážeme. Nech X1, . . . , Xn, n ∈ N je náhodný výber z Gumbelovho rozdelenia s parametrami µ a σ, t.j. jeho distribučná funkcia má tvar G0,µ,σ = exp −e− x−µ σ . Kapitola 2. Metóda blokových maxím 17 Potom distribučná funkcia maxima je podľa (2.1) tvaru Gn(x) = P(Mn ≤ x) = [F(x)]n = exp −e−x−µ σ n Postupne upravujeme Gn(x) = exp −n · e−x−µ σ = exp −eln n e−x−µ σ = exp −e σ ln n−(x−µ) σ = = exp −e− x−(µ+σ ln n) σ = G0,µ∗,σ(x), čo je distribučná funkcia Gumbelovho rozdelenia s parametrom polohy µ∗ = µ+σ ln n posunutým vzhľadom na pôvodné rozdelenie a s rovnakým parametrom merítka σ. Podobne, pre Fréchetovo rozdelenie s distribučnou funkciou s parametrami µ a σ definovanou ako G1,µ,σ(x) = exp − x − µ σ −α , α > 0 počítajme Gn(x) = [F(x)]n = exp − x − µ σ −α n = exp −n x − µ σ −α = = exp − x − µ n1/ασ −α = G1,µ,σ∗ (x). Úpravami sme prišli opäť k distribučnej funkcii Fréchetovho rozdelenia s parametrami µ a σ∗ = σ · n1/α . Nakoniec uvažujme Weibullovo rozdelenie s distribučnou funkciou G2,µ,σ(x) = exp − − x − µ σ −α , α < 0, s parametrami µ a σ. Platí Gn(x) = [F(x)]n = exp − − x − µ σ −α n = exp −n − x − µ σ −α = = exp − − x − µ n1/ασ −α = G2,µ,σ∗ (x). Pritom G2,µ,σ∗ (x) je znova distribučnou funkciou Weibullovho rozdelenia, tentokrát s parametrami µ a σ∗ = n1/α σ. Z výpočtov plynie, že ak má náhodný výber nejaké z rozdelení extrémnych hodnôt, tak potom maximum tohto výberu má opäť príslušné rozdelenie extrémnych hodnôt, s prípadne zmenenými parametrami polohy alebo merítka. Tým sme ukázali, že EV rozdelenia sú max–stabilné. Kapitola 2. Metóda blokových maxím 18 2.3 Modelovanie metódou blokových maxím V nasledujúcej časti sa zameriame na praktickejšie úlohy pri analýze blokových ma- xím. Nech X1, . . . , XN , N ∈ N je náhodný výber z nejakého rozdelenia s distribučnou funkciou F(x). Pre účely metódy rozdeľme dáta do m blokov po n pozorovaní. Celkovo teda máme N = n · m pozorovaní. Ak n · m = N, môžeme túto skutočnosť pre dostatočne veľký počet pozorovaní zanedbať, ak sa počet pozorovaní v blokoch nelíši dramaticky. V každom bloku potom nájdime maximum Mi, i = 1, . . . , m. Blokové maximá tvoria náhodný výber, s ktorým budeme ďalej pracovať. 2.3.1 Odhad parametrov GEV rozdelenia Pre dostatočne veľký počet pozorovaní v jednotlivých blokoch by postupnosť blokových maxím mala konvergovať ku GEV rozdeleniu. Budeme potom ďalej predpokladať, že M1, . . . , Mm je náhodný výber, ktorý môžeme aproximovať GEV rozdelením s parametrami γ, µ a σ. Tieto parametre však nepoznáme a musíme ich odhadnúť vhodnou metódou, napríklad metódou maximálnej vierohodnosti, či metódou vážených momentov. Metóda maximálnej vierohodnosti Jedným z najbežnejšie používaných odhadov je maximálne vierohodný odhad. Podstatou tejto metódy je určenie vierohodnostnej funkcie náhodného výberu z nejakého rozdelenia. Následnou maximalizáciou tejto funkcie získame takzvané maximálne vierohodné odhady parametrov rozdelenia. Hustota GEV rozdelenia je tvaru gγ,µ,σ(x) = 1 σ 1 + γ(x − µ) σ − 1 γ −1 · exp − 1 + γ(x − µ) σ − 1 γ , pre 1 + γ(x−µ) σ > 0, γ = 0, γ ∈ R, µ ∈ R a σ > 0. Potom logaritmus hustoty GEV rozdelenia je ln gγ,µ,σ(x) = − ln σ − 1 γ + 1 ln 1 + γ(x − µ) σ − 1 + γ(x − µ) σ − 1 γ . Teraz môžeme zhotoviť logaritmickú vierohodnostnú funkciu l(γ, µ, σ; M1, . . . , Mm) pre náhodný výber blokových maxím M1, . . . , Mm, a teda l(γ, µ, σ; M1, . . . , Mm) = m i=1 ln gγ,µ,σ(Mi) = = −m ln σ − 1 γ + 1 m i=1 ln 1 + γ(Mi − µ) σ − m i=1 1 + γ(Mi − µ) σ − 1 γ . Kapitola 2. Metóda blokových maxím 19 Príslušné maximálne vierohodnostné odhady γ, µ a σ získame ako argument maxima logaritmickej vierohodnostnej funkcie. Výraz máme maximalizovať pre všetky hodnoty γ ∈ R, µ ∈ R, σ > 0 spĺňajúce podmienku 1 + γ(Mi − µ) σ > 0, i = 1, . . . , m. Explicitné riešenie tejto úlohy neexistuje, a preto sa na výpočet odhadov parametrov používajú numerické metódy. Poznamenajme, že množina prípustných hodnôt odhadovaných parametrov závisí na pozorovaniach, čo znamená, že nie sú splnené podmienky regularity. Práve podmienkami regularity sa zaoberal Smith [21]. Dokázal, že pre hodnoty γ > −0,5 sú odhady asymptoticky nestranné, konzistentné a majú asymptoticky normálne rozdelenie. V prípade −1 < γ ≤ −0,5 maximálne vierohodné odhady existujú, no nutne nemajú štandardné asymptotické vlastnosti. Nakoniec, pre γ < −1 takéto odhady zvyčajne vôbec neexistujú. V praxi to však nie je prekážka, nakoľko prípad, kedy parameter γ < −0,5 nastáva len naozaj zriedkavo. Metóda vážených momentov Metóda pravdepodobnostne vážených momentov je zovšeobecnenie klasickej momentovej metódy pre odhad parametrov. Definujme najskôr všeobecné pravdepodobnostne vážené momenty. Definícia 2.3.1. Nech X je náhodná veličina z nejakého rozdelenia s distribučnou funkciou F(x). Potom veličiny Mp,r,s = E {Xp [F(X)]r [1 − F(X)]s } , sa nazývajú pravdepodobnostne vážené momenty, pričom p, r a s sú reálne čísla. Ak je možné kvantilovú funkciu F−1 (u) = inf {x ∈ R; F(x) ≥ u} k distribučnej funkcii F vyjadriť explicitne, potom môžeme vážené momenty ekvivalente zapísať v tvare Mp,r,s = 1 0 F−1 p Fr (1 − F)s dF. (2.13) Tento tvar býva často vhodnejší pre výpočet vážených momentov. Práve preto ho využijeme na odhad parametrov GEV rozdelenia, pričom Hosking a kol. [12] odporúčajú uvažovať vážené momenty typu βr = M1,r,0 = E {X [F(X)]r } , r = 0, 1, 2. V prípade GEV rozdelenia s distribučnou funkciou Gγ,µ,σ(x) = exp − 1 + γ(x − µ) σ − 1 γ , 1 + γ(x − µ) σ > 0, γ = 0, Kapitola 2. Metóda blokových maxím 20 majú momenty βr tvar βr = 1 r + 1 µ − σ γ [1 − (r + 1)γ Γ(1 − γ)] , γ < 1, γ = 0. (2.14) Vzťah (2.14) odvodíme. Podľa (2.13) platí βr = 1 0 G−1 Gr dG Kvantilová funkcia GEV rozdelenia je tvaru G−1 = µ + σ γ (− ln G)−γ − 1 , a preto pokračujeme βr = 1 0 µ + σ γ (− ln G)−γ − 1 Gr dG. Zavedením substitúcie t = − ln G dostávame βr = ∞ 0 µ + σ γ t−γ − 1 e−t(r+1) dt = = µ − σ γ ∞ 0 e−t(r+1) dt + σ γ ∞ 0 t−γ e−t(r+1) dt = = µ − σ γ (r + 1)−1 + σ γ ∞ 0 t−γ e−t(r+1) dt. Druhá substitúcia u = (r + 1)t nás pomaly dovedie k výsledku: βr = µ − σ γ (r + 1)−1 + σ γ (r + 1)γ−1 ∞ 0 u−γ e−u du, kde integrál vo výraze je gama funkcia2 Γ(1−µ), a preto vzťah (2.14) platí, čím sme odvodili pravdepodobnostne vážené momenty zovšeobecneného rozdelenia extrémnych hodnôt. Rovnice na získanie odhadov γ, µ, σ parametrov γ, µ, σ dostaneme pomocou vážených momentov GEV rozdelenia, pričom pre r = 0 podľa (2.14) platí β0 = µ − σ γ [1 − Γ(1 − γ)] . (2.15) Ďalej pre r = 1 priebežným využitím vyššie odvodeného vzťahu pre β0 dostávame 2 Gama funkcia je definovaná ako Γ(x) = ∞ 0 ux−1 e−u du pre x > 0. Kapitola 2. Metóda blokových maxím 21 β1 = 1 2 µ − σ γ [1 − 2γ Γ(1 − γ)] , 2β1 = µ − σ γ [1 − Γ(1 − γ) + Γ(1 − γ) − 2γ Γ(1 − γ)] , 2β1 = µ − σ γ [1 − Γ(1 − γ)] + σ γ Γ(1 − γ)(2γ − 1), 2β1 − β0 = σ γ Γ(1 − γ)(2γ − 1). Obdobne získame postupnými úpravami rovnicu pre r = 2, a teda β2 = 1 3 µ − σ γ [1 − 3γ Γ(1 − γ)] , 3β2 = µ − σ γ [1 − Γ(1 − γ)] + σ γ Γ(1 − γ)(3γ − 1), 3β2 − β0 = σ γ Γ(1 − γ)(3γ − 1), 3β2 − β0 = 2β1 − β0 2γ − 1 , 3β2 − β0 2β1 − β0 = 3γ − 1 2γ − 1 . Celkovo sme obdržali tri rovnice o troch neznámych β0 = µ − σ γ [1 − Γ(1 − γ)] , 2β1 − β0 = σ γ Γ(1 − γ)(2γ − 1), (2.16) 3β2 − β0 2β1 − β0 = 3γ − 1 2γ − 1 , ktoré je potreba vyriešiť. Pritom βr, r = 0, 1, 2 nahradíme nejakým vhodným odhadom. Napríklad Landwehr [16] odvodil nestranný odhad pre zoradený náhodný výber M(1) ≤ . . . ≤ M(n) v tvare βr = 1 n n j=1 r l=1 j − l n − l M(j), r = 1, 2, pričom β0 = n−1 n j=1 M(j). Potom odhad γ parametru γ získame z poslednej rovnice sústavy (2.16) numerickým riešením. Ďalej, po dosadení tohto odhadu do druhej rovnice, máme σ = γ(2β1 − β0) Γ(1 − γ)(2γ − 1) . Kapitola 2. Metóda blokových maxím 22 Nakoniec, z prvej rovnice získame odhad parametru µ µ = β0 + σ γ [1 − Γ(1 − γ)] . Odhady γ, µ, σ rovnako ako v prípade maximálne vierohodných odhadov nespĺňajú podmienky regularity, avšak sú asymptoticky normálne pre γ < 0,5. Naviac, odhady metódou vážených momentov vykazujú lepšie správanie ako maximálne vierohodné odhady v prípade malých výberov, pozri [12]. Existuje mnoho ďalších typov odhadov pre parametre EV a GEV rozdelení. Často sa využíva Hillov odhad ([2], [7], [17]), ďalej Pickandsov či Dekkers–Einmahl–de Haanov odhad ([7]) alebo klasická metóda momentov ([3]), prípadne metóda vážených najmenších štvorcov ([2]). 2.3.2 Spracovanie dát V nasledujúcej časti sa zameriame na vybrané postupy pre štatistickú analýzu extrémnych hodnôt. Postupne rozoberieme niekoľko otázok – aká je pravdepodobnosť, že náhodný jav nadobudne nejakej extrémnej hodnoty a aká je úroveň a doba návratu extrémneho javu. Pravdepodobnosť nastania extrémnej udalosti Uvažujme príklad, kedy máme k dispozícii dáta o množstve napadnutých zrážok na určitom území. Aká je pravdepodobnosť, že nasledujúci deň napadne napr. 100 mm zrážok? Matematicky formulujme úlohu nasledovne. Nech X1, . . . , XN je náhodný výber z rozdelenia s distribučnou funkciou F(x), pričom pozorovania sú rozdelené do m blokov s počtom pozorovaní n v každom bloku. Odhadnime teraz pravdepodobnosť P(Xi > x) pre nejaké x veľké. Ak je teda dĺžka bloku n pozorovaní, tak podľa (2.1) platí, že distribučná funkcia blokového maxima Mn je P(Mn ≤ x) = [P(Xi ≤ x)]n = [1 − P(Xi > x)]n , pričom v poslednom vyjadrení sa nachádza člen, ktorý je predmetom nášho záujmu a ktorý chceme v určitej rozumnej forme vyjadriť. Zároveň poznáme tvar distribučnej funkcie Gγ,µ,σ(x) zovšeobecneného rozdelenia extrémnych hodnôt (2.3), pomocou ktorého modelujeme blokové maximá, a preto platí Gγ,µ,σ(x) = [1 − P(Xi > x)]n , z čoho vyplýva, že P(Xi > x) = 1 − (Gγ,µ,σ(x)) 1 n . Po dosadení distribučnej funkcie GEV rozdelenia do predchádzajúceho vzťahu, pričom γ = 0 a 1 + γ(x − µ)/σ > 0, výraz upravíme do tvaru Kapitola 2. Metóda blokových maxím 23 P(Xi > x) = 1 − exp − 1 + γ(x − µ) σ − 1 γ 1 n = 1 − exp − 1 n 1 + γ(x − µ) σ − 1 γ . Pre veľký rozsah výberu môžeme výraz aproximovať pomocou prvých dvoch členov Taylorovho polynómu pre rozvoj funkcie ex , kde približne platí, že ex = 1 + x. Parametre GEV rozdelenia nahradíme príslušnými odhadmi. Celkovo teda môžeme vysoké (extrémne) pravdepodobnosti odhadnúť vzťahom P(Xi > x) = 1 n 1 + γ(x − µ) σ − 1 γ , (2.17) kde µ, σ a γ sú vhodné odhady príslušných parametrov. Odhad extrémnych kvantilov – úroveň a doba návratu Opäť na príklade dát o množstve napadnutých zrážok uvažujme situáciu, kedy nás zaujíma, aké extrémne množstvo zrážok sa priemerne vyskytne napr. raz za 10 rokov. Tomuto pojmu sa hovorí úroveň návratu (return level), pričom teda chceme určiť veľkosť extrémneho javu, ktorý bude v priemere nastávať so známou zvolenou frekvenciou. V našom kontexte sa jedná o odhad vysokých (extrémnych) kvantilov GEV rozdelenia. Hľadajme úroveň, ktorá bude i-tym blokovým maximom Mi, i = 1, . . . , m, prekročená raz za k časových okamihov, čo sa dá pomocou definície (1 − 1/k)-kvantilu vyjadriť ako P Mi > q1− 1 k = 1 k , k ∈ N, (2.18) pričom q1−1/k je hľadaná úroveň. Pre väčšiu prehľadnosť ďalšieho výpočtu označme q := q1−1/k. Pre dostatočný rozsah náhodného výberu aproximujme výraz (2.18) pomocou distribučnej funkcie GEV rozdelenia (2.3): 1 − Gγ,µ,σ(q) = 1 − exp − 1 + γ(q − µ) σ − 1 γ = 1 k , z ktorého chceme vyjadriť úroveň návratu q, pričom γ = 0 a 1 + γ(x − µ)/σ > 0. Preto pokračujme vo výpočte Kapitola 2. Metóda blokových maxím 24 exp − 1 + γ(q − µ) σ − 1 γ = 1 − 1 k , − 1 + γ(q − µ) σ − 1 γ = ln 1 − 1 k , 1 + γ(q − µ) σ = − ln 1 − 1 k −γ , odkiaľ q = µ + σ γ − ln 1 − 1 k −γ − 1 . Po nahradení µ, σ a γ prislúchajúcimi odhadmi µ, σ, γ získame q = q1− 1 k = µ + σ γ − ln 1 − 1 k −γ − 1 , (2.19) čo je tzv. odhad úrovne návratu. S extrémnymi kvantilmi súvisí aj ďalší pojem – doba návratu (return period) extrémnej udalosti. V tomto prípade je situácia opačná, snažíme sa napr. určiť, priemerne za aké obdobie denné zrážky prekročia 150 mm. To znamená, že k dispozícii máme údaj o úrovni návratu q = q1−1/k, no priemerná frekvencia, s ktorou bude jav nastávať, nie je známa a snažíme sa ju určitým spôsobom odhadnúť. Zo vzťahu (2.18) P Mi > q1− 1 k = 1 k preto vyjadríme dobu návratu k k = 1 P Mi > q1− 1 k = 1 1 − exp − 1 + γ(q−µ) σ −1/γ Prislúchajúci odhad je potom tvaru k = 1 1 − exp − 1 + γ(q−µ) σ −1/γ . Vzťah medzi dobou a úrovňou návratu graficky znázorňuje tzv. return-level plot. Tvar grafu závisí na hodnote parametru tvaru γ zovšeobecneného rozdelenia extrémnych hodnôt. Ak γ = 0, graf je lineárny. Ak γ > 0, graf má konvexný tvar a ak γ < 0, závislosť je konkávneho tvaru. Na obr. 2.4 môžeme vidieť príklad takéhoto grafu, pričom na vodorovnej osi je znázornená doba návratu a na zvislej úroveň návratu. Kapitola 2. Metóda blokových maxím 25 q q q q q q q q q qqqqqqqqqq qqqqqq qqqqq q q q q q q q q q q q q q q q q q q q q 0.2 0.5 1.0 2.0 5.0 10.0 20.0 50.0 100.0 0246810 Return Level Plot Return Period ReturnLevel Obr. 2.4: Return–level plot pre γ = 0,5. Nakoniec tejto kapitoly zostáva otázkou, ako je vhodné zvoliť počet pozorovaní v jednotlivých blokoch, čo patrí medzi najkomplikovanejšie úlohy v metóde blokových maxím. Tu proti sebe vystupujú dve skutočnosti. Pokiaľ určíme bloky, ktoré obsahujú malý počet pozorovaní, získame tak mnoho blokových maxím a to má za následok zlé asymptotické vlastnosti modelu a jeho veľké vychýlenie. No na druhú stranu, veľký počet pozorovaní v blokoch vedie k malému počtu blokových maxím, čo vytvára model s veľkým rozptylom. V praxi sa doposiaľ nezaužíval žiadny jednotný postup na určenie veľkosti blokov a počtu pozorovaní v nich. Niektorí autori navrhujú využiť Hillov odhad parametru γ, ktorý sa vynesie do grafu v závislosti na počte pozorovaní v jednotlivých blokoch, a ďalej odhadom určiť bod, od ktorého sa graf stabilizuje (pozri napr. [22]). Nevýhodou tohto postupu je akási subjektívnosť každého jedinca pri určovaní bodu stability. Poznamenajme však, že bloky sa v praxi mnohokrát volia tak, aby korešpondovali s časovým obdobím, ku ktorému sa analyzované dáta viažu. Blokové maximá sú potom často štvrťročné, ročné, desaťročné maximá a podobne. Kapitola 3 Peaks-over-threshold V nasledujúcej metóde uvažujme postupnosť i.i.d. náhodných veličín X1, . . . , Xn s distribučnou funkciou F. Viac, ako veľkosť a správanie sa týchto veličín nás však budú zaujímať náhodné veličiny Y1, . . . , YNu , Nu ≤ n, ktoré udávajú veľkosť prekročenia nejakej vopred danej medze u. Prekročenia sú definované ako Yi = Xi − u, Xi ≥ u, pričom Nu udáva počet prekročení hranice u náhodnými veličinami X1, . . . , Xn, matematicky zapísané ako Nu = card {i : i = 1, . . . , n; Xi ≥ u } . Náhodné veličiny Yi, niekedy nazývané aj excesy, má preto zmysel uvažovať, len pokiaľ k prekročeniu danej medze došlo. Zostrojme teraz distribučnú funkciu excesov Fu(x) pre náhodné veličiny Yi, i = 1, . . . , Nu. Zrejme sa preto jedná o podmienenú distribučnú funkciu Fu(x) = P(Yi ≤ x|Xi ≥ u) = P(Xi − u ≤ x|Xi ≥ u). (3.1) Využitím vzťahu pre podmienenú pravdepodobnosť získavame Fu(x) = P(u ≤ Xi ≤ u + x) P(Xi ≥ u) = F(u + x) − F(u) 1 − F(u) , x ≥ 0. (3.2) Môžeme si všimnúť, že podmienená distribučná funkcia Fu(x) závisí na výške zvolenej medze u. V predchádzajúcom odvodení sme uvažovali prekročenia nad určitú medzu, t.j. pre maximálne extrémne udalosti. Analogicky môžeme uvažovať aj extrémne udalosti opačným smerom, čiže udalosti, ktoré (extrémne) klesli pod určitú hranicu. Ďalej sa v celej práci budeme zaoberať situáciami, kedy pozorovanie prekročí vopred danú hranicu smerom nahor. Podobne ako pri metóde blokových maxím v predchádzajúcej kapitole, vhodne normované excesy budú konvergovať k určitej triede rozdelení, čo ukážeme neskôr. Uveďme preto teraz tieto rozdelenia spolu s ich distribučnými funkciami. Definícia 3.0.1. Pre modelovanie excesov používame rozdelenia definované pomocou nasledujúcich distribučných funkcií: – 26 – Kapitola 3. Peaks-over-threshold 27 −1.0 −0.6 −0.2 0.0 0.00.51.01.5 Beta rozdelenie x hustota α 1.5 3 4.5 0 1 2 3 4 5 0.00.51.01.5 Exponenciálne rozdelenie x hustota 1 2 3 4 5 0.00.51.01.5 Paretovo rozdelenie x hustota α 1.5 3 4.5 Obr. 3.1: Hustoty beta, exponenciálneho a Paretovho rozdelenia pre rôzne hodnoty parametru α. Exponenciálne: W0(x) = 1 − e−x , x ≥ 0, Paretovo: W1(x) = 1 − x−α , x ≥ 1, Beta: W2(x) = 1 − (−x)α , −1 ≤ x ≤ 0. kde α > 0 sa nazýva parameter tvaru. Hustoty rozdelení z definície 3.0.1 môžeme vidieť na obrázku 3.1. Exponenciálne a Paretovo rozdelenie sa používajú pri modelovaní maxím a beta rozdelenie pri modelovaní miním. Všeobecný tvar rozdelení dostaneme zavedením parametru polohy µ a parametru merítka σ, t.j. Wi,µ,σ(x) = Wi x − µ σ , µ ∈ R, σ > 0, i = 0, 1, 2. Rozdelenia z definície 3.0.1 môžeme definovať aj jednotným kompaktnejším spôsobom, pomocou zovšeobecneného Paretovho rozdelenia. Definícia 3.0.2. Zovšeobecnené Paretovo (GP - Generalised Pareto) rozdelenie má distribučnú funkciu v tvare Wγ(x) = 1 − (1 + γx)− 1 γ , γ ∈ R, γ = 0, pričom γ sa nazýva parameter tvaru. Hustota zovšeobecneného Paretovho rozdelenia pre rôzne hodnoty parametru γ je zobrazená na obr. 3.2. Nakoľko GP rozdelenie má ľavý koncový bod rovný nule, Kapitola 3. Peaks-over-threshold 28 0.0 0.5 1.0 1.5 2.0 2.5 0.00.51.01.52.02.5 GPD x hustota γ −1.2 −0.8 0 0.8 Obr. 3.2: Hustota zovšeobecneného Paretovho rozdelenia s rôznymi hodnotami parametru γ. parameter polohy µ v našom prípade nepotrebujeme uvažovať. Úplnú triedu GP rozdelenia potom získame pridaním parametru merítka σ, t.j. Wγ,σ(x) = 1 − 1 + γx σ − 1 γ , γ ∈ R, σ > 0. Poznamenajme, že GP a GEV rozdelenia majú medzi sebou úzky vzťah, ktorý sa dá jednoducho vyjadriť ako Wγ(x) = 1 + ln Gγ(x). Ďalej sa môžeme presvedčiť, že pre rôzne hodnoty parametru tvaru γ je GP rozdelenie len pozmeneným zápisom rozdelení z definície 3.0.1. Pre γ = 0 sa jedná o exponenciálne rozdelenie, distribučnú funkciu GP rozdelenia však musíme spojito dodefinovať v nule, a teda previesť limitný prechod. Preto lim γ→0 Wγ(x) = lim γ→0 1 − (1 + γx)− 1 γ = 1 − lim γ→0 (1 + γx)− 1 γ = 1 − e−x , čo je distribučná funkcia exponenciálneho rozdelenia. V prípade γ > 0 a γ = 1/α, α > 0 pre distribučnú funkciu GP rozdelenia získavame Wγ(x) = 1 − (1 + γx)− 1 γ = 1 − 1 + x α −α = W1 1 + x α , čím sme prišli k distribučnej funkcii Paretovho rozdelenia v bode 1 + x/α. Nakoniec, pre γ < 0 a zároveň γ = −1/α, α > 0 máme Wγ(x) = 1 − (1 + γx)− 1 γ = 1 − 1 − x α α = W2 −1 + x α , Kapitola 3. Peaks-over-threshold 29 čiže sme získali beta rozdelenie v bode −1+x/α, a tým aj overili korektnosť tvrdenia o zovšeobecnenom Paretovom rozdelení. Zaoberajme sa teraz limitným rozdelením excesov. Špeciálne nás zaujíma, aké rozdelenie budú mať so zvyšujúcou sa hranicou u. Preto teraz uvedieme akúsi analógiu Fisherovej-Tippetovej vety 2.0.2 pre excesy. Veta 3.0.1 (Balkemova-de Haanova, Pickandsova). Nech X1, . . . , Xn je náhodný výber z nejakého rozdelenia. Ak existujú postupnosti konštánt an > 0, bn ∈ R tak, že Fu(anx + bn) u ↑xF −→ H(x), pre nejakú spojitú distribučnú funkciu H(x), potom platí lim u↑xF |Fu(anx + bn) − Wγ,σ(u)(x)|= 0, x ∈ R, kde Wγ,σ(u)(x) je distribučná funkcia GP rozdelenia. Dôkaz. Dôkaz je opäť náročný, je inšpirovaný Gnedenkovým dôkazom Fisherovej-Tippetovej vety a ako prví ho nezávisle na sebe sformulovali Balkema a de Haan ([1]) a Pickands ([18]). Prechádzajúca veta nám vraví, že jediné možné limitné rozdelenie excesov náhodného výberu pre u ↑ xF 1 je zovšeobecnené Paretovo rozdelenie. Budeme teda predpokladať, že excesy Yi = Xi −u (za podmienky Xi ≥ u) majú asymptoticky GP rozdelenie Yi ≈ Wγ,σ a budeme ich modelovať pomocou tohto rozdelenia, pričom parametre γ a σ odhadneme z dát. Ukázať konvergenciu na konkrétnom príklade, ako sme to urobili v predchádzajúcej kapitole, už nie je také jednoduché, nakoľko sa v limite má zvyšovať nielen počet pozorovaní náhodného výberu, ale aj výška thresholdu u. Tak isto z predchádzajúcej vety vidíme, že variabilita σ(u) závisí na voľbe hranice u. Čím nižšie threshold zvolíme, tým bude variabilita väčšia a opačne. S tým súvisí problém, ako vhodne zvoliť hranicu u, ktorému sa budeme venovať nižšie. Medzi GP a GEV rozdelením je, ako už bolo povedané, blízky vzťah. Dokonca platí, že rozdelenia, pre ktoré normované maximá konvergujú k GEV rozdeleniu, tvoria množinu, pre ktorú príslušné rozdelenie excesov konverguje s rastúcou medzou k zovšeobecnenému Paretovmu rozdeleniu. Zaujímavé pritom je, že parameter tvaru γ nie je závislý na voľbe thresholdu u a je v oboch prípadoch rovnaký (pozri [7]). 3.1 Stabilita rozdelení pre modelovanie excesov V predchádzajúcej kapitole sme ukázali, že EV rozdelenia sú max-stabilné. Podobnú vlastnosť majú aj rozdelenia excesov, je však nutné pozmeniť definíciu max-stability, nakoľko nás nezaujíma stabilita blokových maxím, ale stabilita prekročení nad hranicu u. 1 u ↑ xF značí, že u sa blíži k xF zdola. Kapitola 3. Peaks-over-threshold 30 Definícia 3.1.1. Rozdelenie náhodného výberu X1, . . . , Xn sa nazýva stabilné voči excesom, ak pre jeho nedegenerovanú distribučnú funkciu F(x) platí, že existuje vhodná postupnosť konštánt an > 0 a bn ∈ R takých, že Fu(anx + bn) = F(x) pre všetky x ∈ R, pričom Fu(x) distribučná funkcia excesov definovaná podľa (3.1). Túto vlastnosť overíme najskôr na samotnom GP rozdelení. Nech preto náhodný výber X1, . . . , Xn je zo zovšeobecneného Paretovho rozdelenia s distribučnou funkciou Wγ,σ(x) = 1 − 1 + γx σ − 1 γ . Pre pevne zvolené u a excesy Yi = Xi −u (za podmienky Xi ≥ u), i = 1, . . . , n podľa (3.2) platí Fu(x) = P(u ≤ Xi ≤ u + x) P(Xi ≥ u) = Wγ,σ(u + x) − Wγ,σ(u) 1 − Wγ,σ(u) . Dosadením a malou úpravou získavame Fu(x) = 1 + γ σ u − 1 γ − 1 + γ σ (u + x) − 1 γ 1 + γ σ u − 1 γ = 1 − 1 + γ σ (u + x) 1 + γ σ u − 1 γ , pričom úpravou na spoločného menovateľa dostaneme Fu(x) = 1 − σ + γ(u + x) σ + γu − 1 γ = 1 − 1 + γx σ + γx − 1 γ = Wγ,σ∗ (x), čo je opäť distribučná funkcia GP rozdelenia s rovnakým parametrom tvaru γ a zmeneným parametrom σ∗ = σ + γu. Tým sme ukázali, že excesy majú znova GP roz- delenie. Už sme ukázali, že GP rozdelenie je len iným zápisom rozdelení na modelovanie excesov z definície (3.0.1), čo automaticky znamená, že tieto rozdelenia sú taktiež stabilné voči excesom. Môže byť však zaujímavé pozrieť sa, ako sa zmenia jednotlivé parametre polohy a merítka v týchto konkrétnych prípadoch. Prejdime preto postupne ku všetkým trom rozdeleniam. Nech náhodný výber X1, . . . , Xn je z exponenciálneho rozdelenia s distribučnou funkciou W0,µ,σ(x) = 1 − exp − x − µ σ , µ ∈ R, σ > 0. Potom podmienená distribučná funkcia excesov má tvar Kapitola 3. Peaks-over-threshold 31 Fu = W0,µ,σ(u + x) − W0,µ,σ(u) 1 − W0,µ,σ(u) = exp −u−µ σ − exp −u+x−µ σ exp −u−µ σ , a teda postupnými úpravami prídeme až k Fu(x) = 1 − exp −u+x−µ σ exp −u−µ σ = 1 − exp − u + x − µ σ + u − µ σ = = 1 − exp − x σ = W0,µ∗,σ. Tým sme prišli k distribučnej funkcii exponenciálneho rozdelenia s parametrom polohy µ∗ = 0 a rovnakým parametrom merítka σ. Táto vlastnosť exponenciálneho rozdelenia je známa aj v iných oblastiach a aplikáciach matematiky, kedy hovoríme, že exponenciálne rozdelenie je tzv. bez pamäti. Podobne, pre náhodný výber z Paretovho rozdelenia s distribučnou funkciou W1,µ,σ(x) = 1 − x − µ σ −α , α > 0 môžeme odvodiť Fu(x) = W1,µ,σ(u + x) − W1,µ,σ(u) 1 − W1,µ,σ(u) = u−µ σ −α − u +x−µ σ −α u−µ σ −α , čo po úpravách vedie ku konečnému výsledku Fu(x) = 1 − x − (µ − u) u − µ −α = W1,µ∗,σ∗ (x). Toto ukazuje, že excesy nad určitú medzu z náhodného výberu z Paretovho rozdelenia majú tak isto Paretovo rozdelenie, no so zmenou parametrov, kde parameter polohy je µ∗ = µ − u a parameter merítka σ∗ = u − µ. Nakoniec, táto vlastnosť platí aj pre náhodný výber z beta rozdelenia s distribučnou funkciou W2,µ,σ = 1 − − x − µ σ α , α > 0. Obdobným postupom získavame Fu(x) = W2,µ,σ(u + x) − W2,µ,σ(u) 1 − W2,µ,σ(u) = −u−µ σ α − −u+x−µ σ α −u−µ σ α . Úpravou tohto výrazu sa získame Fu(x) = 1 − − x − (µ − u) µ + u α = W2,µ∗,σ∗ (x). Tým overujeme platnosť stability voči excesom aj pre beta rozdelenie. V tomto prípade sa zmenil parameter polohy na µ∗ = µ − u a parameter merítka na σ∗ = µ + u. Všetky rozdelenia sú preto stabilné voči excesom. Kapitola 3. Peaks-over-threshold 32 3.2 Metóda POT pri praktickom spracovaní dát Zopakujme, že uvažujeme náhodný výber X1, . . . , Xn s distribučnou funkciou F. Pri použití POT metódy pri modelovaní dát je našou úlohou najskôr odhadnúť parametre GP rozdelenia, ďalej si vhodne zvoliť prah u a definovať excesy. Potom môžeme skúmať rôzne vlastnosti v závislosti na požadovaných výsledkoch, čomu sú venované nasledujúce časti práce. 3.2.1 Odhad parametrov GP rozdelenia Podľa vety 3.0.1 vieme, že pre dostatočne veľké u môžeme excesy Yi = Xi−u, Xi ≥ u, i = 1, . . . , Nu aproximovať zovšeobecneným Paretovým rozdelením. Predpokladajme teda, že náhodný výber Y1, . . . , YNu je náhodný výber z GP rozdelenia s distribučnou funkciou Wγ,σ(x). Parametre γ a σ však nepoznáme, a teda ich musíme odhadnúť z dát. V nasledujúcej časti sú popísané dve metódy zhodné s metódami použitými na odhad parametrov GEV rozdelenia vyššie, a to metóda maximálnej vierohodnosti a metóda vážených momentov. Metóda maximálnej vierohodnosti Podstata tejto metódy je obdobná s prípadom z predchádzajúcej kapitoly, postupujme preto analogicky. Hustotu GP rozdelenia získame deriváciou jeho distribučnej funkcie, a preto wγ,σ(x) = 1 σ 1 + γx σ − 1 γ −1 , γ ∈ R, σ > 0. Logaritmus hustoty teda vyzerá nasledovne ln wγ,σ(x) = − ln σ − 1 γ + 1 ln 1 + γx σ . Potom logaritmická vierohodnostná funkcia bude v tvare l(γ, σ; Y1, . . . , YNu ) = −Nu ln σ − 1 γ + 1 Nu i=1 ln 1 + γYi σ , (3.3) pričom tento výraz maximalizujeme pre σ > 0, γ ∈ R a za podmienky 1 + γYi σ > 0, i = 1, . . . , Nu. Maximalizáciu je potrebné previesť priamo numericky, a príslušné maximálne vierohodné odhady γ a σ získame ako argument maxima logaritmickej vierohodnostnej funkcie (3.3). Poznamenajme, že pri odhadovaní parametrov nie sú splnené podmienky regularity. Avšak pre hodnoty γ > −0,5 táto metóda funguje klasicky a odhady sú asymptoticky nestranné, konzistentné a majú asymptoticky normálne rozdelenie (pozri [8]). Kapitola 3. Peaks-over-threshold 33 Metóda vážených momentov Odhad parametrov GP rozdelenia pomocou metódy pravdepodobnostne vážených momentov prevedieme obdobne ako pri GEV rozdelení. Pripomeňme, že pravdepodobnostne vážené momenty sú definované ako Mp,r,s = E {Xp [F(X)]r [1 − F(X)]s } , kde F je distribučná funkcia náhodnej veličiny X. V prípade odhadu parametrov zovšeobecneného Paretovho rozdelenia odporúčajú Hosking a Wallis (pozri [13]) uvažovať vážené momenty v tvare αs = M1,0,s = E {X [1 − F(X)]s } , s = 0, 1. Toto môžeme podľa (2.13) vyjadriť ako αs = 1 0 F−1 (1 − F)s dF, (3.4) kde F−1 je kvantilová funkcia k distribučnej funkcii F. Pripomeňme, že GP rozdelenie má distribučnú funkciu Wγ,σ(x) = 1 − 1 + γx σ − 1 γ , γ ∈ R, σ > 0. Potom kvantilová funkcia bude W−1 = σ γ (1 − W)−γ − 1 . Dosadením tejto funkcie do (3.4) získame αs = 1 0 W−1 (1 − W)s dW = 1 0 σ γ (1 − W)−γ − 1 (1 − W)s dW V integráli zavedieme substitúciu t = 1 − W, čím získame riešenie nasledovne αs = σ γ 1 0 t−γ − 1 ts dt = σ γ 1 s − γ + 1 − 1 s + 1 a teda αs = σ (s − γ + 1)(s + 1) . Nakoľko chceme odhadnúť dva parametre, potrebujeme vyjadriť α0 a α1, čiže α0 = σ 1 − γ , α1 = σ 2(2 − γ) . Kapitola 3. Peaks-over-threshold 34 Úpravou týchto rovníc získame odhady γ a σ parametrov γ a σ v tvaroch γ = 4α1 − α0 2α1 − α0 , σ = 2α0α1 α0 − 2α1 , pričom αs sa nahradí nejakým odhadom. Napríklad Landwehr navrhol v [16] nestranný odhad pre zoradený náhodný výber Y(1), . . . , Y(Nu) v tvare αs = 1 n Nu j=1 s k=1 Nu − j − k + 1 Nu − k Y(j). Tieto odhady opäť nespĺňajú podmienky regularity. Avšak pre hodnotu parametru γ > −0,5 je v [13] ukázané, že v takomto prípade majú odhady štandardné vlastnosti nestrannosti, konzistentnosti a majú taktiež asymptoticky normálne rozdelenie. 3.2.2 Spracovanie dát Analogicky ako pri metóde blokových maxím opísanej v predchádzajúcej kapitole, sa teraz zamerajme na vybrané štatistické úlohy, ktoré nás môžu pri analýze dát zaujímať. Odhadneme pravdepodobnosť, že nastane extrémna udalosť a pozrieme sa na úroveň a dobu návratu extrémnej udalosti v prípade uvažovania POT metódy. Pravdepodobnosť nastania extrémnej udalosti Predstavme si, že je našou úlohou odhadnúť pravdepodobnosť, že si klient poisťovne bude nárokovať nejaké veľmi vysoké poistné plnenie. Túto otázku môžeme riešiť nasledovne. Nech X1, . . . , Xn je náhodný výber z rozdelenia s distribučnou funkciou F. Zaujíma nás pravdepodobnosť P(Xi > x), pre nejaké veľké x. Pravdepodobnosť v prípade POT metódy uvažujeme ako podmienenú a pre prípad, že došlo k prekročeniu medze u, t.j. ak x ≥ u. Preto podľa definície podmienenej pravdepodobnosti máme P(Xi > x|Xi ≥ u) = P(Xi > x) P(Xi ≥ u) . Odtiaľ vyjadríme a postupne upravujeme P(Xi > x) = P(Xi ≥ u) · P(Xi > x|Xi ≥ u) = = P(Xi ≥ u) · (1 − P(Xi ≤ x|Xi ≥ u)) , pričom zavedením excesu Yi = Xi − u, i = 1, . . . , Nu, Xi ≥ u získame Kapitola 3. Peaks-over-threshold 35 P(Xi > x) = P(Xi ≥ u) · (1 − P(Yi ≤ x − u|Xi ≥ u)) = P(Xi ≥ u) · (1 − Fu(x − u)) , kde Fu(x − u) je podmienená distribučná funkcia excesov. Podľa vety 3.0.1 môžeme pre dostatočne veľkú hranicu u aproximovať Fu(x−u) pomocou GP rozdelenia s distribučnou funkciou Wγ,σ(x − u). Preto P(Xi > x) = P(Xi ≥ u) · 1 + γ(x − u) σ − 1 γ . (3.5) Potom parametre γ a σ nahradíme ich príslušnými odhadmi γ a σ získanými nejakou vhodnou metódou. Ďalej odhadnime P(Xi ≤ u) výrazom Nu/n, kde Nu je počet pozorovaní, ktoré prekročili danú medzu u a n je počet všetkých pozorovaní. Celkovo získavame vzťah P(Xi ≥ x) = Nu n 1 + γ(x − u) σ − 1 γ , (3.6) ktorý môžeme použiť na odhad extrémnych pravdepodobností. Odhad extrémnych kvantilov – úroveň a doba návratu Obdobne ako v predchádzajúcej metóde, úroveň návratu je vlastne veľkosť extrémneho javu, ktorý bude priemerne nastávať s vopred zvolenou frekvenciou. Zmena však nastáva v chápaní plynutia frekvencie. Kým pri metóde blokových maxím sme ju merali v blokoch (t.j. typicky to je ročná, desaťročná frekvencia a pod.), pri POT metóde sa frekvencia meria v pozorovaniach, čo už častokrát nie je také intuitívne. V nasledujúcom texte budeme uvažovať úroveň návratu, ktorá bude prekročená extrémnym javom raz za k pozorovaní. Poznamenajme však, že ak by nás zaujímala úroveň návratu za K rokov, môžeme zaviesť substitúciu k = K · ny, kde ny je počet pozorovaní za rok. Toto môžeme využiť v prípade, že ny je známe a v jednotlivých rokoch sa príliš nelíši (pri veľkom počte pozorovaní môžeme malé rozdiely zanedbať). Hľadajme teraz úroveň návratu q := q1−1/k, ktorá bude i-tym pozorovaním prekročená raz za k pozorovaní ako (1 − 1/k)-kvantil, t.j. podľa definície P(Xi > q) = 1 k , k ∈ N. Podľa vzťahu (3.5), ktorý sme odvodili vyššie máme P(Xi ≥ u) · 1 + γ(q − u) σ − 1 γ = 1 k , Kapitola 3. Peaks-over-threshold 36 z čoho úpravami vyjadríme požadovaný kvantil q = u + σ γ [(k · P(Xi ≥ u))γ − 1] . Potom príslušný odhad bude q = q1− 1 k = u + σ γ k Nu n γ − 1 , (3.7) kde γ, σ sú vhodné odhady parametrov GP rozdelenia γ, σ a Nu/n je odhad pre P(Xi ≥ u), pričom Nu je počet pozorovaní, ktoré prekročili medzu u a n je celkový počet pozorovaní. Ďalšou charakteristikou, ktorá nás zaujíma, je doba návratu extrémnej udalosti. Je to vlastne priemerná frekvencia výskytu extrémneho javu pri známej úrovni návratu q = q1−1/k, t.j. snažíme sa nejakým spôsobom odhadnúť k. Vychádzajme zo vzťahu pre (1 − 1/k)-kvantil rozdelenia P(Xi > q) = 1 k . Potom k = 1 P(Xi > q) (3.5) = 1 P(Xi ≥ u) 1 + γ(q−u) σ . Príslušný odhad k je teda tvaru k = Nu n 1 + γ(q − u) σ −1 , (3.8) pričom γ, σ sú vhodné odhady γ, σ, Nu je počet prekročení nad medzu u a n je celkový počet pozorovaní. Voľba thresholdu u Voľba thresholdu u je jednou z najkritickejších úloh POT analýzy. Ak zvolíme u príliš vysoko, získame model s málo excesmi a tým pádom s veľkým rozptylom. Naopak, pri nízko zvolenej hranici u vzniká vychýlený model. Je preto potrebné nájsť nejaký kompromis. K nájdeniu hranice u sa používajú viaceré metódy. Jednou z nich je grafická metóda založená na linearite strednej hodnoty excesov pre GP rozdelenie, ktorú si najskôr definujme. Definícia 3.2.1. Stredná hodnota excesov (mean excess function) je funkcia definovaná ako e(u) = E(X − u|X > u), u ∈ R, kde u je zvolená hranica prekročenia. Kapitola 3. Peaks-over-threshold 37 Poznamenajme, že ak má náhodná veličina Z zovšeobecnené Paretovo rozdelenie s parametrami γ a σ, potom stredná hodnota tejto veličiny je EZ = σ 1 − γ γ < 1. Ďalej sme v kapitole 3.1 odvodili, že excesy Yi = Xi −u za podmienky Xi ≥ u, majú zovšeobecnené Paretovo rozdelenie s parametrami γ a σ∗ = σ + γu. Skombinovaním týchto dvoch informácií dostaneme pre strednú hodnotu excesov GP rozdelenia vzťah e(u) = σ + γu 1 − γ , čo je lineárna funkcia premennej u. Toto je charakteristická vlastnosť GP rozdelenia, v prípade uvažovania iného rozdelenia by sme získali zložitejšie vyjadrenie. V praxi, pre náhodný výber X1, . . . , Xn sa funkcia e(u) môže odhadnúť výberovým priemerom excesov, t.j. e(u) = 1 Nu Nu i=1 Yi, kde Yi je exces, ak k nemu došlo a Nu je počet prekročení hranice u. Následne vynesieme do grafu dvojice bodov X(i), e(X(i)) , i = 1, . . . , n , čím vznikne tzv. mean excess plot a X(1), . . . , X(n) je pritom usporiadaný náhodný výber. Pokiaľ pozorovania pochádzajú z GP rozdelenia, potom by sa mal graf po určitom čase stabilizovať a správať sa približne lineárne. Potom hodnotu thresholdu u vyberieme na základe toho, od ktorého pozorovania sa graf javí lineárny. Nevýhodou je subjektivita jednotlivca pri výbere hranice touto optickou metódou, nakoľko interpretácia tohto grafu je častokrát náročná. Ďalšou prekážkou môže byť fakt, že stabilizácia grafu v skutočnosti častokrát vôbec nenastane, nakoľko sa rozdelenie excesov nepodarí dostatočne aproximovať GP rozdelením, čo robí úlohu určenia hranice u ešte komplikovanejšou. Niektoré ďalšie metódy určenia thresholdu u boli navrhnuté, pozri [20]. Jednou z možností je určiť hranicu prekročenia u ako nejaký vysoký kvantil, často sa používa 90 % alebo 96 % kvantil. Ďalšou alternatívou je označiť ako u hodnotu pozorovania X(n−j+1), zoradeného náhodného výberu X(1) ≤ . . . ≤ X(n), pričom j = √ n alebo j = n2/3 ln(ln(n)). Pri praktickom spracovaní dát je vhodné zvoliť viacero hodnôt pre threshold u a výsledné modely vzájomne porovnať. Kapitola 4 Analýza dát Analýza extrémnych udalostí sa uplatňuje v mnohých oblastiach. V nasledujúcej kapitole ju aplikujeme na dáta z oblasti poisťovníctva. Našim cieľom bude na základe známych historických dát poskytnúť odhad extrémnej udalosti a iné vyššie zmieňované charakteristiky. Odpovede na tieto otázky sú v tejto oblasti dôležité, nakoľko pomáhajú poistným spoločnostiam vykonať rozhodnutie napríklad o alokácii kapitálu. Na záver budú výsledky rôznych metód porovnané. Dáta boli modelované v jazyku R , pričom zdrojové kódy je možné nájsť v elektronickom archíve tejto práce. 4.1 Dáta Dáta sú tvorené informáciami o poistných nárokoch americkej poistnej spoločnosti z oblasti súkromného poistenia automobilov [9]. Premenná, ktorá je ďalej v práci využitá obsahuje údaje o veľkosti čiastky, ktorá bola spoločnosťou poistenému vyplatená. Priebeh tejto premennej obsahujúcej 6 773 pozorovaní je zobrazený na obr. 4.1. 0 1000 2000 3000 4000 5000 6000 7000 0200004000060000 Index pozorovania Výškapoistnýchnárokov Obr. 4.1: Výška poistných nárokov. – 38 – Kapitola 4. Analýza dát 39 Minimálna vyplatená čiastka je 9,5 $, maximálna 60 000 $. Ďalšie základné charakteristiky sú zhrnuté v tabuľke 4.1. Priemer 1 853 Rozptyl 7 006 129 Medián 1 002 Smerodajná odchýlka 2 649 1. kvartil 524 Šikmosť 6,23 3. kvartil 2 137 Špicatosť 84,25 Tabuľka 4.1: Základné charakteristiky dátového súboru. Šikmosť je kladná, preto očakávame dáta zošikmené doprava a hodnota špicatosti naznačuje ťažké chvosty rozdelenia dát. Tieto vlastnosti vizuálne potvrdzuje aj histogram dát na obr. 4.2 (s hodnotami vyššími ako 25 000 odseknutými). Histogram poistných nárokov Výska poistných nárokov Frekvencia 0 5000 10000 15000 20000 25000 050010001500 Obr. 4.2: Histogram výšky poistných nárokov. Vidíme, že dataset je výrazne zošikmený doprava s ťažkými chvostami a tým pádom by modelovanie pomocou normálneho rozdelenia nebolo vhodné. Využijeme preto metódy na modelovanie extrémnych situácií popísané v predchádzajúcich kapitolách práce. 4.2 Modelovanie metódou blokových maxím Kľúčovým začiatočným rozhodnutím je výber počtu blokov, ktorý bude pre náš model použitý. Vzhľadom na to, že neexistuje nejaká zaužívaná metóda pre túto voľbu, prevedieme ju porovnaním vlastností modelov s rôznymi počtami blokov. Predovšetkým sa upriamime na to, ako dobre je výber blokových maxím aproximovaný GEV rozdelením na základe grafických nástrojov. Následne sa môžeme zamerať na odhad vybraných charakteristík. Kapitola 4. Analýza dát 40 m = 20 Výška poistného nároku Hustota 0 10000 30000 50000 0.000000.000040.000080.00012 MLE PWM m = 65 Výška poistného nároku Hustota 0 10000 30000 500000.000000.000040.000080.00012 MLE PWM m = 100 Výška poistného nároku Hustota 0 10000 30000 50000 0.000000.000040.000080.00012 MLE PWM Obr. 4.3: Histogram blokových maxím pre počet blokov m = 20, 65 a 100. Červená krivka je hustota odhadnutá na základe maximálne vierohodných odhadov a modrá krivka je odhadnutá hustota na základe pravdepodobnostne vážených momentov. Výber počtu blokov Modely boli porovnávané okrem iného podľa histogramov pre rôzne počty blokových maxím. Aj na základe tejto metódy sme dospeli k rozhodnutiu, že najvhodnejšie bude zvoliť počet blokov m = 65, čím sa vytvárajú bloky s veľkosťou n = 104 (resp. 105) pozorovaní. Pre ilustráciu môžeme na obr. 4.3 pozorovať histogramy pre dáta rozdelené do 20, 65 a 100 blokov. Zároveň sú v grafoch znázornené odhadnuté hustoty zovšeobecneného rozdelenia extrémnych hodnôt - červenou krivkou vzhľadom na parametre odhadnuté maximálnou vierohodnosťou (MLE) a modrou pravdepodobnostne váženými momentami (PWM). Tak isto boli porovnané p-hodnoty testu zhody empirického rozdelenia blokových maxím s GEV rozdelením s neznámymi parametrami. Využitá bola technika založená na výberovom korelačnom koeficiente, ktorá bola navrhnutá pre rozdelenia stabilné voči maximám (kam náš prípad spadá) v publikácii [11]. Výsledky sú zhrnuté v tabuľke 4.2. Podľa tohto testu na hladine významnosti 5 % nezamietame nulovú hypotézu, že rozdelenie maxím pre počet blokov 20 a 65 pochádza z GEV rozdelenia. Naopak, pre 100 blokových maxím túto hypotézu zamietame. m p-hodnota 20 0,156 65 0,188 100 0,029 Tabuľka 4.2: P-hodnoty prislúchajúce testu zhody s GEV rozdelením pre rôzny počet blokov m. Kapitola 4. Analýza dát 41 Histogram Výška poistného nároku Hustota 0 10000 30000 50000 0.000000.000040.00008 10000 20000 30000 40000 100002000030000400005000060000 Q−Q plot Teoretické kvantilyEmpirickékvantily 0.0 0.2 0.4 0.6 0.8 1.0 0.00.20.40.60.81.0 P−P plot Teoretická distribucná funkcia Empirickádistribucnáfunkcia Obr. 4.4: Histogram, Q-Q plot a P-P plot pre maximálne vierohodné odhady parametrov GEV rozdelenia. Odhad parametrov GEV rozdelenia Odhady parametrov pre výber blokových maxím Mi, i = 1, . . . , 65 boli prevedené metódou maximálnej vierohodnosti a metódou pravdepodobnostne vážených momentov. Pre nami vybraný počet blokov sú ich hodnoty zhrnuté v tabuľke 4.3, pričom sa oba prípady od seba vzájomne markantne nelíšia. γ µ σ MLE 0,30 12 269 4 023 PWM 0,33 12 239 3 919 Tabuľka 4.3: Odhady parametrov GEV rozdelenia metódou maximálnej vierohodnosti (MLE) a metódou pravdepodobnostne vážených momentov (PWM). Pozrime sa aj na grafické metódy pre posúdenie vhodnosti vybraného modelu, ktoré sú založené na porovnávaní teoretického a empirického rozdelenia. Q-Q plot porovnáva teoretické a empirické kvantily rozdelenia a P-P plot teoretickú a empirickú distribučnú funkciu. Pretože sú tieto diagnostické grafy približne lineárne, môžeme povedať, že GEV rozdelenie aproximuje naše dáta relatívne dobre. Na obr. 4.4 vidíme tieto grafy pre maximálne vierohodné odhady a na obr. 4.5 pre pravdepodobnostne vážené momentové odhady. V oboch prípadoch je približná linearita relatívne dobre zachovaná, menšiu komplikáciu ale predstavujú dve pozorovania, ktoré sú odľahlejšie ako ostatné, čo je badateľné napr. z Q-Q plotu. To môže budúce odhady pravdepodobností do istej miery skresliť. Extrémne pravdepodobnosti V predchádzajúcej časti práce sme odvodili, aký je odhad pravdepodobnosti, že nasledujúce pozorovanie bude mať vyššiu hodnotu ako nejaké veľké x, čo odhadujeme pomocou vzťahu (2.17). V našom prípade to znamená pravdepodobnosť, že nasledujúci zákazník, ktorý do poisťovne príde, bude požadovať poistný nárok vyšší ako Kapitola 4. Analýza dát 42 Histogram Výška poistného nároku Hustota 10000 30000 50000 0.000000.000040.00008 10000 20000 30000 40000 100002000030000400005000060000 Q−Q plot Teoretické kvantilyEmpirickékvantily 0.0 0.2 0.4 0.6 0.8 1.0 0.00.20.40.60.81.0 P−P plot Teoretická distribucná funkcia Empirickádistribucnáfunkcia Obr. 4.5: Histogram, Q-Q plot a P-P plot pre parametrov GEV rozdelenia pomocou pravdepodobnostne vážených momentov. x dolárov. Tieto pravdepodobnosti sú v prípade skúmaných dát relatívne nízke pre vysoké čiastky poistného plnenia, napr. pravdepodobnosť, že nasledujúci zákazník si bude nárokovať viac ako 35 000 $ odhadujeme približne na 0,0349 %. Tieto charakteristiky teda nepatria medzi najužitočnejšie. Preto je napríklad zaujímavejšie pozorovať, aká je pravdepodobnosť, že medzi nasledujúcimi c zákazníkmi príde aspoň jeden zákazník, ktorý požaduje poistný nárok vyšší ako x $. To odhadneme nasledovnou úvahou: pravdepodobnosť, že aspoň jeden zákazník bude mať nárok vyšší ako x je doplnok k pravdepodobnosti, že žiadny z c zákazníkov nebude mať nárok vyšší ako x, čo sa dá vyjadriť ako 1 − P(X1 ≤ x) · . . . · P(Xc ≤ x) = 1 − [1 − P(Xi > x)]c , nakoľko v našom modeli predpokladáme, že X1, . . . , Xc je postupnosť i.i.d. náhodných veličín. Pravdepodobnosť P(Xi > x) už potom odhadneme známym vzťahom (2.17) odvodeným v teoretickej časti práce. Na obrázku 4.6 vidíme odhady týchto pravdepodobností pre celú škálu výšky poistných nárokov a rôzne počty uvažovaných novoprichádzajúcich klientov ci, ktorí majú tieto poistné nároky. Dve krivky reprezentujú výsledok výpočtu za použitia maximálne vierohodných odhadov GEV rozdelenia a odhadov získaných pomocou pravdepodobnostne vážených momentov. V prvom prípade pre c1 = 1 000 klientov sú pritom takmer na nerozoznanie a relatívne presne sa kopírujú. V ostatných príkladoch pre počty klientov c2 = 5 000 a c3 = 10 000 rozdiely mierne rastú. Na grafoch taktiež pozorujeme, že čím vyšší počet klientov c uvažujeme, tým je krivka reprezentujúca odhad pravdepodobnosti menej „strmá a viac posunutá doprava, teda odhadovaná pravdepodobnosť sa s vyšším počtom klientov taktiež zvyšuje. Je potrebné poznamenať, že pre nízke poistné nároky treba brať dané grafy s rezervou, je dobré však vidieť aspoň približný priebeh. Metódy modelovania extrémnych udalostí sú navrhnuté pre vysoké x, t.j. vysoké hodnoty požadovaných poistných nárokov a tým pádom v opačnom prípade dobre nefungujú. Na otázku, čo už sú extrémne poistné plnenia nie je jednoznačná odpoveď, no vzhľadom na povahu dát by mohli byť Kapitola 4. Analýza dát 43 odhady relevantné približne od čiastky 40 000 dolárov. 0 50000 100000 150000 0.00.20.40.60.81.0 c1 = 1 000 Výška poistného nároku Pravdepodobnost MLE PWM 0 50000 100000 150000 0.00.20.40.60.81.0 c2 = 5 000 Výška poistného nároku Pravdepodobnost MLE PWM 0 50000 100000 150000 0.20.40.60.81.0 c3 = 10 000 Výška poistného nároku Pravdepodobnost MLE PWM Obr. 4.6: Odhad extrémnych pravdepodobností - pravdepodobnosť, že medzi nasledujúcimi c klientmi bude aspoň jeden s poistným nárokom vyšším ako x dolárov. MLE - maximálne vierohodné odhady, PWM - odhady metódou pravdepodobnostne vážených momentov pre parametre GEV rozdelenia. V tabuľke 4.4 vidíme konkrétny príklad odhadnutých pravdepodobností. Označme pritom pci , i = 1, 2, 3 ako odhad pravdepodobnosti, že medzi nasledujúcimi ci zákazníkmi sa objaví aspoň jeden s nárokom vyšším ako x $. Počty zákazníkov ci sú pritom opäť c1 = 1 000, c2 = 5 000 a c3 = 10 000. Aj číselne sa potvrdili zvyšujúce sa rozdiely odhadov pravdepodobností použitím MLE a PWM metódy pre rastúci počet uvažovaných prichádzajúcich zákazníkov. Interpretácia napríklad druhého riadku tabuľky (a) je, že pravdepodobnosť, že medzi nasledujúcimi 1 000 zákazníkmi sa objaví aspoň jeden klient s nárokom vyšším ako 60 000 $ je 6 %, že sa objaví medzi nasledujúcimi 5 000 zákazníkmi je 26 % a medzi nasledujúcimi 10 000 zákazníkmi je 45 %, všetko pri použití maximálne vierohodných odhadov GEV rozdelenia. Teda potvrdilo sa zrejmé, čím vyšší počet prichádzajúcich zákazníkov budeme uvažovať, tým vyššia je pravdepodobnosť, že sa medzi nimi objaví aspoň nejaký s požiadavkou na vysoké poistné plnenie. x pc1 pc2 pc3 40 000 0,20 0,68 0,90 60 000 0,06 0,26 0,45 80 000 0,02 0,11 0,21 100 000 0,01 0,05 0,10 (a) x pc1 pc2 pc3 40 000 0,22 0,72 0,92 60 000 0,07 0,31 0,52 80 000 0,03 0,14 0,26 100 000 0,02 0,08 0,14 (b) Tabuľka 4.4: Odhad extrémnych pravdepodobností (a) maximálne vierohodnými odhadmi (b) odhadmi pomocou pravdepodobnostne vážených momentov GEV rozde- lenia. Kapitola 4. Analýza dát 44 Úroveň a doba návratu V prípade dát o poistných nárokoch je zaujímavé preskúmať, aká veľkosť poistného nároku bude prekonaná v priemere jedenkrát za c zákazníkov. Inak povedané, hľadáme úroveň návratu, ktorá bude priemerne nastávať so zvolenou frekvenciou. Matematicky teda odhadujeme extrémny kvantil GEV rozdelenia podľa vzťahu (2.19). Jediná odlišnosť tkvie v tom, že kým sme v predchádzajúcom prípade počítali s periódou k ekvivalentnou dĺžke jedného bloku, tentokrát sa zameriame na frekvenciu v podobe počtu zákazníkov. Preto k vyjadríme ako k = c n , kde c je počet zákazníkov a n je počet pozorovaní v jednom bloku. V tabuľke 4.5 sú odhadnuté úrovne návratu q pre rôzny počet zákazníkov pre obe metódy odhadov parametrov rozdelenia, ktoré sú popísané vyššie. Vidíme napríklad, že priemerne raz za 10 000 klientov sa objaví taký, ktorý si bude od poisťovne nárokovať poistné plnenie vyššie ako 51 239, resp. 54 094 dolárov. Ako je zrejmé, úrovne návratu s rastúcim počtom zákazníkov taktiež stúpajú. Odhad úrovne návratu pomocou pravdepodobnostne vážených momentov je pritom vyšší ako maximálne vierohodný, pričom ich rozdiel navzájom sa prehlbuje. Kým pre 5 000 zákazníkov sa odhadnuté úrovne líšia o približne 4 %, pre 75 000 zákazníkov sú odlišné o 10 %. c qMLE qPWM 5 000 41 396 42 989 10 000 51 239 54 094 20 000 63 304 68 038 35 000 75 016 81 871 50 000 83 560 92 122 75 000 94 436 105 345 Tabuľka 4.5: Odhad úrovne návratu pre rôzny počet zákazníkov c. Možnosťou takisto je, že poisťovňa bude naopak viac sústredená na priemernú frekvenciu nejakej hodnoty poistného plnenia. V tomto prípade sa preto zaujímame o dobu návratu extrémnej udalosti, opäť vyjadrenú pomocou počtu zákazníkov. Doby návratu pre rôzne výšky poistného plnenia q a pre obe metódy odhadov parametrov GEV rozdelenia sú vypísané v tabuľke 4.6. Môžeme napríklad pozorovať, že poistné plnenie vyššie ako 60 000 $ sa vyskytne priemerne raz za 43 292, resp. 32 637 zákaz- níkov. V odhadoch oboch charakteristík – úrovne aj doby návratu – môžeme sledovať, že GEV rozdelenie odhadnuté metódou pravdepodobnostne vážených momentov priraďuje extrémnym udalostiam vyššiu váhu v porovnaní s maximálne vierohodne odhadnutým GEV rozdelením, pričom rozdiely sa stávajú pomerne vysoké s rastúcou hodnotou požadovaného poistného nároku. Inak povedané, prvé z menovaných má ťažšie chvosty voči druhému. Tento vzťah môžeme vidieť aj na obr. 4.7, pričom tento graf sa nazýva return-level plot. Zobrazuje závislosť medzi úrovňou a dobou návratu, ktorá je uvažovaná pomocou zákazníkov. Kapitola 4. Analýza dát 45 q cMLE cPWM 40 000 4 476 4 025 60 000 16 770 13 676 80 000 43 292 32 637 100 000 90 690 64 078 Tabuľka 4.6: Odhad doby návratu pre rôzne veľkosti poistného plnenia q.02000060000100000 Return−level plot Doba návratu Úrovennávratu 200 500 1000 2000 5000 10000 20000 50000 MLE PWM MLE PWM Obr. 4.7: Return-level plot, kde MLE značí maximálne vierohodné odhady a PWM odhady pomocou pravdepodobnostne vážených momentov GEV rozdelenia. 4.3 Modelovanie POT metódou Druhým prístupom je modelovanie dát pomocou metódy prekročenia nad určitú medzu. Podobne, ako bolo v predchádzajúcom prípade kritickou úlohou vybrať vhodný počet blokov, do ktorých pozorovania rozdelíme, teraz je našim problémom nájsť vhodnú medzu. Po vyriešení tohto problému budú odhadnuté parametre GP rozdelenia a prevedená podobná analýza ako v prípade modelovania pomocou blokových maxím. Voľba thresholdu V kapitole 3 boli navrhnuté viaceré prístupy, ako sa k dopracovať k vhodnej voľbe hranice prekročenia u. Jedným z nich bolo vybrať takú hodnotu, od ktorej sa graf – tzv. mean excess plot – stabilizuje. Tento graf znázorňujúci závislosť strednej hodnoty excesov na výške thresholdu vidíme na obr. 4.8. Nejaví však očividné známky stabilizácie, ako sa teoreticky predpokladá. Napriek tomu sme zvolili dve veľkosti tejto medze, ktoré boli ďalej preskúmavané, pričom vybrané hodnoty sa nachádzajú na horizontálnej osi a sú znázornené prerušovanou čiarou. Výška thresholdu odpovedajúca zvoleným hodnotám je 14 300, reps. 18 500 dolárov. Takisto skúmané boli hranice zvolené na základe spomenutých heuristických odhadov. Modely boli opäť vyhodnocované najmä na základne grafických metód. Porovnanie histogramov s vykreslenou hustotou GP rozdelenia v troch prípadoch môžeme Kapitola 4. Analýza dát 46 0 10000 20000 30000 40000 50000 60000 050001500025000 Mean excess plot u e^(u) Obr. 4.8: Mean excess plot. pozorovať na obr. 4.9. Vidíme, že MLE a PWM odhady sú takmer identické a od seba nerozoznateľné vo všetkých vybraných možnostiach. V tabuľke 4.7 sa potom nachádzajú počty excesov prislúchajúce výškam hranice u, ktoré sú na obrázku zobrazené. Vybrané thresholdy odpovedajú (zľava) hodnote zvolenej na základe mean excess plotu, 96 % kvantilu a pozorovaniu X(n−j+1), kde j = n2/3 ln(ln(n)), pričom predpokladáme, že vstupné dáta tvoria náhodný výber X1, . . . , Xn. u 14 300 7210 3875 Nu 47 271 780 Tabuľka 4.7: Počet excesov Nu pre rôzne výšky hranice prekročenia u. Kombináciou všetkých faktorov bola zvolená výška thresholdu ako u = 7 210, čím získavame 271 excesov. Nakoniec bol prevedený test zhody empirického rozdelenia dát a zovšeobecneného Paretovho rozdelenia. Vybraný test bol navrhnutý v [23], pričom testovanie hypotézy prebieha metódou bootstrap. Výsledná p-hodnota je rovná 0,25, a preto na hladine významnosti 5 % nezamietame hypotézu, že výber excesov pochádza z GP rozdelenia. Odhad parametrov GP rozdelenia Predpokladajme, že prekročenia Yi, i = 1, . . . , Nu, získané voľbou hranice u = 7 210 tvoria náhodný výber zo zovšeobecneného Paretovho rozdelenia. Odhady jeho parametrov γ a σ metódou maximálnej vierohodnosti a pravdepodobnostne váženými momentami sú zhrnuté v tabuľke 4.8. Obe metódy produkujú veľmi blízke odhady, dokonca odhady parametru γ sa v tomto prípade líšia až na vzdialenejších desatinných pozíciách. Vizuálne metódy posúdenia vhodnosti modelu pre maximálne vierohodné odhady GEV rozdelenia sú zobrazené na obrázku 4.10. Vidíme, že GP rozdelenie relatívne Kapitola 4. Analýza dát 47 u = 14300 Exces Hustota 0 10000 30000 0.000000.000100.000200.00030 MLE PWM u = 7210 Exces Hustota 0 10000 30000 500000.000000.000100.000200.00030 MLE PWM u = 3875 Exces Hustota 0 10000 30000 50000 0.000000.000100.000200.00030 MLE PWM Obr. 4.9: Histogram excesov v porovnaní s hustotou GP rozdelenia. Parametre GP rozdelenia boli odhadnuté maximálne vierohodnostnou (MLE) metódou a pravdepodobnostne váženými momentami (PWM). γ σ MLE 0,279 2 960 PWM 0,281 2 950 Tabuľka 4.8: Odhad parametrov GP rozdelenia pomocou metódy maximálnej vierohodnosti a pomocou pravdepodobnostne vážených momentov. dobre aproximuje výber excesov, slabinou môžu byť dve odľahlé pozorovanie, ktoré sú zreteľné na Q-Q plote. Z dôvodu veľkej podobnosti odhadov parametrov pri použití oboch metód sú obdobné diagnostické grafy pre metódu pravdepodobnostne vážených momentov vynechané. Extrémne pravdepodobnosti Ak by poisťovňu zaujímala pravdepodobnosť, že nasledujúci zákazník, ktorý bude mať voči nej poistný nárok, bude požadovať čiastku vyššiu ako x dolárov, mohla by ju odhadnúť pomocou vzťahu (3.6). Konkrétnejšie, pravdepodobnosť, že nasledujúci klient bude požadovať poistné plnenie vyššie ako 35 000 $ by sa podľa neho odhadla na 0,0398%. Pri metóde blokových maxím sme však odhadovali pravdepodobnosti, že medzi nasledujúcimi prichádzajúcimi c zákazníkmi sa objaví aspoň jeden klient s poistným nárokom vyšším ako nejaké vysoké x. Analogicky riešme túto úlohu metódou POT. Obdobne, túto pravdepodobnosť vyjadríme ako 1 − (1 − P(Xi > x))c kde P(Xi > x) sa odhadne pomocou vzťahu (3.6). Vývoj tejto veličiny v závislosti na výške extrémneho poistného nároku je zobrazený na obr. 4.11, pričom uvažujeme Kapitola 4. Analýza dát 48 Histogram Výška excesu Hustota 0 10000 30000 50000 0.000000.000100.000200.00030 0 10000 20000 30000 40000 0100003000050000 Q−Q plot Teoretické kvantilyEmpirickékvantily 0.0 0.2 0.4 0.6 0.8 1.0 0.00.20.40.60.81.0 P−P plot Teoretická distribucná funkcia Empirickádistribucnáfunkcia Obr. 4.10: Histogram, Q-Q plot a P-P pre maximálne vierohodné odhady parametrov GP rozdelenia. rôzny počet prichádzajúcich zákazníkov ci, i = 1, 2, 3. Zároveň znova upozorňujeme, že odhady pravdepodobností pre malé poistné čiastky sú skôr orientačné, nakoľko metóda na tieto situácie nie je stavaná, pretože uvažuje len extrémne udalosti. Tabuľka 4.9 potom ukazuje niektoré vybrané odhady pravdepodobností pre rôzne poistné nároky. Napr. pravdepodobnosť, že medzi nasledujúcimi c2 = 5 000 zákazníkmi príde aspoň jeden zákazník s nárokom vyšším než 60 000 dolárov je odhadnutá ako 0,28 pre maximálne vierohodné odhady GP rozdelenia. Odhady pravdepodobností za použitia MLE a PWM metódy sa líšia úplne minimálne. x pc1 pc2 pc3 40 000 0,23 0,72 0,92 60 000 0,06 0,28 0,48 80 000 0,02 0,12 0,22 100 000 0,01 0,06 0,11 (a) x pc1 pc2 pc3 40 000 0,23 0,73 0,93 60 000 0,07 0,29 0,49 80 000 0,03 0,12 0,22 100 000 0,01 0,06 0,11 (b) Tabuľka 4.9: Odhad extrémnych pravdepodobností (a) maximálne vierohodnými odhadmi (b) odhadmi metódou pravdepodobnostne vážených momentov GP rozdelenia. Úroveň a doba návratu Zaoberajme sa teraz opäť rovnakými úlohami ako pri predchádzajúcej metóde blokových maxím. Hľadajme preto najskôr úroveň návratu poistného nároku, ktorý bude v priemere nastávať s danou frekvenciou. Výhodou oproti prvej metóde je fakt, že kým pri blokových maximách bola chápaná frekvencia v rámci blokov a bolo treba zaviesť substitúciu, aby sme mohli úroveň návratu vyjadriť pomocou zákazníkov, teraz nič podobné nie je potrebné. Frekvencia je chápaná rovno v pozorovaniach, t.j. v zákazníkoch, a môžeme preto na odhad priamo využiť vzťah (3.7). Výsledné Kapitola 4. Analýza dát 49 0 50000 100000 150000 0.00.20.40.60.81.0 c1 = 1 000 Výška poistného nároku Pravdepodobnost MLE PWM 0 50000 100000 150000 0.00.20.40.60.81.0 c2 = 5 000 Výška poistného nároku Pravdepodobnost MLE PWM 0 50000 100000 150000 0.00.20.40.60.81.0 c3 = 10 000 Výška poistného nároku Pravdepodobnost MLE PWM Obr. 4.11: Odhad extrémnych pravdepodobností - pravdepodobnosť, že medzi nasledujúcimi c klientmi bude aspoň jeden s poistným nárokom vyšším ako x dolárov. MLE - maximálne vierohodné odhady, PWM - pravdepodobnostne vážené odhady GP rozdelenia. odhady úrovne návratu q pre rôzny počet zákazníkov c sú viditeľné v tabuľke 4.10. Rozdiely medzi použitím MLE a PWM odhadov nie sú markantné. c qMLE qPWM 5 000 43 133 43 294 10 000 53 062 53 325 20 000 65 110 65 518 35 000 76 688 77 253 50 000 85 068 85 756 75 000 95 665 96 519 Tabuľka 4.10: Odhad úrovne návratu pre rôzny počet zákazníkov c. Druhým pohľadom je snaha odhadnúť dobu návratu k. Využijeme na to opäť vzťah (3.8) odvodený v predchádzajúcich kapitolách. Pýtame sa, priemerne raz za koľko zákazníkov sa vyskytne poistné plnenie vyššie ako q dolárov? Súhrnné odhady pozorujeme v tabuľke 4.11. To znamená, že do poisťovne by mal prísť zákazník s poistným nárokom vyšším ako napr. 80 000 dolárov v priemere jedenkrát za približne 40 000 zákazníkov. Vzťah medzi úrovňou a dobou návratu opäť pekne ilustruje return-level plot na obr. 4.12. Rozdiel v týchto odhadoch použitím MLE a PWM odhadov GP rozdelenia je nepatrný. Kapitola 4. Analýza dát 50 q cMLE cPWM 40 000 3 895 3 853 60 000 15 149 14 860 80 000 40 471 39 431 100 000 87 442 84 718 Tabuľka 4.11: Odhad doby návratu pre rôzne veľkosti poistného plnenia q.2000060000100000 Return−level plot Doba návratu Úrovennávratu 100 200 500 1000 2000 5000 10000 20000 50000 MLE PWM MLE PWM Obr. 4.12: Return-level plot, kde MLE značí maximálne vierohodné odhady a PWM odhady GP rozdelenia pomocou pravdepodobnostne vážených momentov. 4.4 Porovnanie metódy blokových maxím a POT metódy Pre lepšiu prehľadnosť sú v ďalšom texte zhrnuté výsledky z analýzy dát metódou blokových maxím aj POT metódou. Pozrime sa najskôr na samotný odhad parametru tvaru pre obe metódy. Podľa kapitoly 3 by za splnenia všetkých teoretických predpokladov mal byť odhad parametru γ rovnaký pre GEV aj GP rozdelenie. Jeho odhady vo všetkých prípadoch, t.j. maximálne vierohodný odhad aj odhad pomocou pravdepodobnostne vážených momentov, môžeme vidieť v tabuľke 4.12. MLE PWM GEV 0,297 0,332 GP 0,279 0,281 Tabuľka 4.12: Odhad parametru γ. Poznamenajme, že odhady nie sú príliš rozdielne, avšak tieto výsledky by mali byť brané orientačne, nakoľko práve splnenie teoretických predpokladov uvažujeme, no nemôžeme si byť istý, že sa tak aj naozaj udialo. Kapitola 4. Analýza dát 51 Extrémne pravdepodobnosti Čo sa odhadu extrémnych pravdepodobností týka, obe metódy vykazujú podobné výsledky. Ako príklad uveďme odhad pravdepodobnosti, že do poisťovne príde klient s poistným nárokom vyšším ako 35 000 dolárov. Pri metóde blokových maxím ju odhadujeme na 0,0349 % a pri využití POT metódy na 0,0398%, takže vidíme, že sa líšia až na treťom desatinnom mieste. Ďalej sme skúmali pravdepodobnosť, že medzi nasledujúcimi c zákazníkmi sa objaví aspoň jeden s poistným nárokom vyšším ako x dolárov. Výsledky týchto odhadov pri použití maximálne vierohodných odhadov rozdelení sa pre obe metódy líšia len mierne, ako môžeme pozorovať na obrázku 4.13. Rozdiely sa zrejme prehlbujú tým, čím vyšší počet prichádzajúcich klientov c uvažujeme. Pritom metóda blokových maxím je v tomto prípade o niečo konzervatívnejšia ako POT metóda a priraďuje vysokým čiastkam o čosi nižšie pravdepodobnosti. 0 40000 80000 120000 0.00.20.40.60.81.0 c1 = 1 000 Výška poistného nároku Pravdepodobnost BM POT 0 40000 80000 120000 0.00.20.40.60.81.0 c2 = 5 000 Výška poistného nároku Pravdepodobnost BM POT 0 40000 80000 120000 0.00.20.40.60.81.0 c3 = 10 000 Výška poistného nároku Pravdepodobnost BM POT Obr. 4.13: Odhad extrémnych pravdepodobností, BM – metódou blokových maxím, POT – metódou peaks-over-threshold. Lepšiu ilustráciu poskytuje obrázok 4.14, na ktorom sú zobrazené odhady extrémnych pravdepodobností oboma spôsobmi modelovania a taktiež pre obe metódy odhadu parametrov rozdelení. Tu vidíme, že najvyššie odhady pravdepodobností získavame použitím metódy blokových maxím a odhadov GEV rozdelenia získaných pomocou pravdepodobnostne vážených momentov. Naopak, ako najnižie sú odhadované pri použití maximálne vierohodných odhadov pri metóde blokových maxím a niekde v strede sú odhady týchto pravdepodobností POT metódou v oboch prípadoch. Pritom posledné dve spomenuté sa asi od čiastky 100 000 dolárov ustálujú na približne rovnakých hodnotách. Úroveň a doba návratu Porovnanie úrovne a doby návratu pre metódy blokových maxím a peaks-over-threshold je najjednoduchšie pomocou return-level plotu, ktorý vidíme na obr. 4.15. Kapitola 4. Analýza dát 52 40000 60000 80000 100000 120000 0.00.20.40.60.81.0 c = 10 000 Výška poistného nároku Pravdepodobnost Blokové maximá MLE PWM Peaks−over−threshold MLE PWM Obr. 4.14: Odhad extrémnych pravdepodobností. Znova sú v ňom vykreslené obe metódy s odhadmi príslušných rozdelení nadobudnutými metódou maximálnej vierohodnosti aj pomocou pravdepodobnostne vážených momentov. Získavame pritom podobné výsledky ako v prechádzajúcom prípade. Najmenšiu priemernú dobu návratu extrémneho poistného plnenia odhaduje model blokových maxím pri použití GEV rozdelenia s odhadmi pomocou pravdepodobnostne vážených momentov, pričom rozdiel oproti ostatným trom spôsobom narastá so zvyšujúcim sa uvažovaným poistným plnením. Napríklad, ako vidíme v tabuľke 4.13, pre konkrétnu hodnotu poistného nároku vo výške 100 000 dolárov získavame odhad priemernej doby návratu o približne 26 000 zákazníkov nižší pri použití PWM odhadov zovšeobecneného rozdelenia extrémnych hodnôt ako pri použití maximálne vierohodných odhadov, čo už vytvára značný rozdiel. Odhady pomocou POT metódy sa opäť držia medzi spomínanými odhadmi metódou blokových maxím. MLE PWM GEV 90 690 64 078 GP 87 442 84 718 Tabuľka 4.13: Odhad doby návratu pre hodnotu poistného nároku q = 100 000. V prípade analyzovaných dát sa ako praktickejšia metóda na modelovanie javí peaks-over-threshold. Metóda blokových maxím môže totižto z analýzy veľmi ľahko vynechať dôležité dáta. Ak sa v jednom bloku aj nachádza viac udalostí s extrémnym dopadom, ďalej uvažované pri modelovaní bude len blokové maximum s najvyššou hodnotou. Naopak, do modelu sa môžu zaniesť blokové maximá, ktoré nie sú až také významné. Vzhľadom na povahu dát o poistných nárokoch, kedy v rámci nich neuvažujeme žiadne časové periódy, by sme preto zvolili metódu POT, ktorá dostupné Kapitola 4. Analýza dát 53 4000080000 Return−level plot Doba návratu Úrovennávratu 2000 5000 10000 20000 40000 70000 Blokové maximá MLE PWM Peaks−over−threshold MLE PWM Obr. 4.15: Return-level plot. dáta využíva lepšie. Naviac, odhad zovšeobecneného Paretovho rozdelenia pri použití MLE a PWM odhadov je v oboch prípadoch veľmi blízky, preto ani nevzniká dilema, ktorý z týchto dvoch odhadov uprednostniť. Avšak použitie metódy blokových maxím môže byť v určitých situáciach taktiež preferované. Či už ide o stav, kedy sú k dispozícii len maximá nejakej premennej (ako napr. ročné maximálne teploty) alebo práve v prípade, kedy sa voľba veľkosti bloku vyskytuje prirodzene (dáta sú mesačné, ročné, desaťročné a pod.) 4.5 Parametrický model Na záver práce sa pokúsme v stručnosti modelovať dáta klasickým parametrickým modelom a odhadnúť pritom pravdepodobnosť, že do poisťovne príde klient s vysokým poistným nárokom ako aj úroveň a dobu návratu vysokého poistného plnenia. Rozdiel oproti vyššie analyzovaným metódam je v tom, že tentokrát sa pokúsime aproximovať celý dátový súbor jedným pravdepodobnostným rozdelením. Analyzovaná dáta sú, ako je uvedené na začiatku tejto kapitoly, silno zošikmené doprava a teda s ťažkými chvostami. Preto sa pri ich ďalšom spracovaní pokúsime využiť niektoré z rozdelení s podobnými vlastnosťami. Do úvahy napríklad pripadajú exponenciálne, Paretovo, lognormálne, či Weibullovo1 . Po preskúmaní a porovnaní možností sme sa rozhodli na modelovanie využiť Weibullovo rozdelenie. Histogram dát s hustotou Weibullovho rozdelenia odhadnutou na základe maximálne vierohodných odhadov vidíme na obrázku 4.16, pričom hodnoty, ktoré presahujú 25 000 dolárov sú odrezané. Ani toto rozdelenie nie je na modelovanie ideálne, no spomedzi ostatných variant sa javilo ako najprijateľnejšia alternatíva. Samotné odhady pravdepodobností objavenia sa zákazníka s vysokým poistným 1 Weibullovym rozdelením v tomto prípade myslíme „klasické Weibullovo rozdelenie, odlišné od toho používaného v kapitole 2. Kapitola 4. Analýza dát 54 Histogram poistných nárokov Výška poistných nárokov Hustota 0 5000 10000 15000 20000 25000 0.00000.00030.0006 Obr. 4.16: Histogram výšky poistných nárokov s odhadnutou hustotou Weibullovho rozdelenia. nárokom sú veľmi nízke až nulové. Napr. pravdepodobnosť, že do poisťovne príde zákazník, ktorý od nej požaduje viac ako 35 000 dolárov je 8,7·10−6 %. Pripomeňme, že pri použití teórie extrémnych hodnôt sme v rovnakom prípade získali odhady rádovo štyrikrát vyššie. Okrem toho sa taktiež pozrime na odhad pravdepodobnosti pc, že medzi nasledujúcimi c = 10 000 zákazníkmi sa objaví aspoň jeden s poistným nárokom vyšším ako x. Tieto pravdepodobnosti tak isto nie sú vysoké a rýchlo klesajú. Pre konkrétne hodnoty sú tieto odhady zhrnuté v tabuľke 4.14. x pc 20 000 0,48453 25 000 0,06814 30 000 0,00770 35 000 0,00087 40 000 0,00010 45 000 0,00001 Tabuľka 4.14: Odhad pravdepodobnosti pc pre c = 10 000 klientov. Dobrú vizuálnu predstavu môžeme získať taktiež z grafu 4.17. Vidíme, že klasický parametrický model priraďuje pravdepodobnostiam, že do poisťovne príde medzi nasledujúcimi c zákazníkmi klient s vysokým poistným nárokom výrazne nižšie odhady v porovnaní s tými, ktoré sme získali pomocou teórie extrémnych hodnôt. Tam, kde ešte len odhady pomocou teórie extrémnych hodnôt ešte len začínajú byť zaujímavé, t.j. okolo hodnoty poistného nároku vo výške približne 40 000, tak sú už odhady pomocou parametrického modelu takmer, či úplne, nulové. Kapitola 4. Analýza dát 55 40000 60000 80000 100000 120000 0.00.20.40.60.81.0 c = 10 000 Výška poistného nároku Pravdepodobnost GP GEV Weibull Obr. 4.17 Úroveň a doba návratu Porovnajme teraz úroveň a dobu návratu extrémnej udalosti medzi parametrickým modelov a výsledkami získanými pomocou teórie extrémnych hodnôt. Úroveň návratu, t.j. veľkosť poistného nároku q, ktorá bude priemerne prekročená raz za c zákazníkov odhadneme ako (1 − 1/c)-kvantil Weibullovho rozdelenia. Pre konkrétne hodnoty a maximálne vierohodné odhady Weibullovho rozdelenia sú výsledky uvedené v tabuľke 4.15. Úroveň návratu stúpa s rastúcim počtom zákazníkov len pomaly. Zatiaľ čo výška poistného nároku, ktorá bude priemerne prekročená raz za 5 000 zákazníkov je odhadnutá na 17 561 dolárov, pre pätnásťkrát viac zákazníkov je to len o zhruba 6 000 dolárov vyššia čiastka. c q 5 000 17 561 10 000 19 089 20 000 20 625 35 000 21 870 50 000 22 666 75 000 23 573 Tabuľka 4.15: Odhad úrovne návratu. Na druhú stranu a v súlade s očakávaniami, odhad doby návratu c rapídne stúpa. Dobu návratu c chápeme ako priemerný počet novoprichádzajúcich zákazníkov, za ktorý sa vyskytne poistný nárok vo výške q. Odhady pre určité zvolené hodnoty sú zhrnuté v tabuľke 4.16. Už len rozdiel medzi odhadom doby návratu pre čiastku 25 000 a 30 000 je celkom priepastný. Vidíme, že tento model prakticky s vysokými poistnými nárokmi nepočíta, resp. uvažuje ich výskyt len s minimálnou pravdepo- Kapitola 4. Analýza dát 56 04000080000 Return−level plot Doba návratu Úrovennávratu 2000 5000 10000 20000 40000 70000 GP GEV Weibull Obr. 4.18: Return-level plot. dobnosťou. Na obr. 4.18 vidíme return-level plot, kde je porovnaný vzťah medzi úrovňou a dobou návratu pre parametrický model a oba modely teórie extrémnych hodnôt. GEV pritom označuje model pomocou zovšeobecneného rozdelenia extrémnych hodnôt a teda metódou blokových maxím a GP značí zovšeobecnené Paretovo rozdelenie, ktoré reprezentuje modelovanie metódou peaks-over-threshold. Z grafu môžeme pozorovať, že krivka zobrazujúca závislosť odhadnutej doby a úrovne návratu parametrickým modelom síce stúpa, no výrazne pomalšie ako ostatné spomenuté. q c 20 000 15 091 25 000 141 702 30 000 1 293 999 40 000 101 087 031 45 000 870 037 866 Tabuľka 4.16: Odhad doby návratu. Parametrický model, ktorý aproximuje vstupné dáta pomocou Weibullovho rozdelenia, výrazne podhodnocuje odhady vysokých pravdepodobností, úrovní a dôb návratu v porovnaní s metódami využívanými v analýze extrémnych udalostí. Analýzou dát sme potvrdili, že modelovaním pomocou metód študovaných v priebehu práce získavame zaujímavé výsledky, ktoré sa z dlhodobého hľadiska môžu javiť ako užitočné. Po zvážení záverov by mohli v našom konkrétnom prípade poslúžiť danej poisťovni napríklad pri posúdení, či rezervy, ktoré si vytvára, sú dostatočné na pokrytie mimoriadne vysokých poistných nárokov, ktoré sú síce vzácne, no s o to väčším dopadom. Záver V tejto diplomovej práci sme sa venovali teórii extrémnych hodnôt a možnostiach jej využitia pri modelovaní dát. V prvých kapitolách sme v ucelenej podobe položili teoretický základ pre túto problematiku. Popísané sú dve široko využívané štatistické metódy pre analýzu extrémnych udalostí. Jedným prístupom je využitie metódy blokových maxím. Podstata metódy tkvie v určení maxím pevne určeného počtu neprekrývajúcich sa blokov a následnom modelovaní tejto série získaných pozorovaní. Predpokladáme pritom, že výber blokových maxím patrí do sféry príťažlivosti zovšeobecneného rozdelenia extrémnych hodnôt, a teda dá sa týmto rozdelením aproximovať. No spomínaný postup je náchylný k neefektívnemu využitiu dostupných informácii, nakoľko z každého bloku sa vyberá len jedna hodnota. Metóda peaks-over-threshold túto nevýhodu odstraňuje. Myšlienkou je zvoliť si hranicu a ďalej modelovať tzv. excesy, t.j. výšky prekročenia nad túto hranicu, ak k tejto udalosti došlo. Potom predpokladáme, že limitné rozdelenie excesov sa dá aproximovať zovšeobecneným Paretovym rozdelením. Metódy boli aplikované na dátach, ktoré sú tvorené informáciami o výške poistných nárokov z automobilového priemyslu. Skúmané boli odhady extrémnych pravdepodobností, úrovne a doby návratu extrémnej udalosti a využité pritom boli maximálne vierohodné a pravdepodobnostne vážené momentové odhady príslušných rozdelení. Vo všetkých prípadoch neboli výsledky analýz markantne rozdielne. Avšak prípad modelovania pomocou blokových maxím s odhadmi GEV rozdelenia pomocou pravdepodobnostne vážených momentov sa zdá byť ako najmenej konzervatívny, pričom priraďuje extrémnym udalostiam najvyššie odhady spomínaných charakteristík spomedzi všetkých možností. Problematickými úlohami pri modelovaní dát boli najmä voľba počtu blokov a výšky thresholdu. V oboch prípadoch bolo toto rozhodnutie zaťažené subjektívnym úsudkom. Porovnaním výsledkov získaných pomocou teórie extrémnych hodnôt a klasického parametrického modelu sme dospeli k záveru, že parametrický model značne podhodnocuje odhady extrémnych pravdepodobností. Je preto dobré uvažovať nad tým, že vezmeme do úvahy práve výsledky analýzy extrémnych situácii, najmä ak by sme sa v danej oblasti chceli pripraviť aj na nečakané udalosti v rámci dlhšieho obdobia. Zaujímavým rozšírením by mohlo byť skúmanie viacrozmerných metód na modelovanie extrémnych udalostí. Príkladom použitia môže byť situácia, kedy okrem výšky klientských poistných nárokov uvažujeme ako druhú zložku ostatné náklady spojené s poistným nárokom, ako napr. náklady spojené s overením pravosti nároku, prípadne právne poplatky. Tento prístup je študovaný napríklad v publikáciach [5] – 57 – Závěr 58 a [6]. Ďalšou možnosťou je skúmanie extrémov stacionárnych časových radov, na čo je poukázané napríklad v [2] alebo [4]. Nakoľko pozorovania v týchto sekvenciách častokrát vykazujú závislosť či sezónnosť, je potrebné základné metódy modifikovať, prípadne vyvinúť ďalšie nástroje pre ich lepšiu analýzu. Tieto a mnohé ďalšie metódy sa v súčasnosti stále rozvíjajú. Seznam použité literatury [1] BALKEMA A.A. a de HAAN L. Residual Life Time of Great Age. The Annals of Probability, 1974, 2(5), pp. 792–804. [2] BEIRLANT, J. a kol. Statistics of Extremes: Theory and Applications. Chichester: John Wiley & Sons, c2004. ISBN 0-471-97647-4. [3] CHRISTOPEIT, N. Estimating parameters of an extreme value distribution by the method of moments. J. Statist. Plann. Inference. 1994, (41), pp. 173–186. [4] COLES, S. An Introduction to Statistical Modeling of Extreme Values. London: Springer, 2001. ISBN 1-85233-459-2. [5] COLES, S. a J. TAWN. Modelling Extreme Multivariate Events. Journal of the Royal Statistical Society. Series B (Methodological). 1991, 53(2), pp. 377–-392. [6] DUPUIS, D. a B. JONES. Multivariate Extreme Value Theory And Its Usefulness In Understanding Risk. North American Actuarial Journal. 2006, 10(4). [7] EMBRECHTS, P., C. KLÜPPELBERG a T. MIKOSCH. Modelling extremal events for insurance and finance. New York: Springer, 1997. ISBN 3540609318. [8] FRANKE, J., W.K. Härdle a C.M. HAFNER. Statistics od Financial Markets. Berlin: Springer, 2011. ISBN 978-3-642-16520-7. [9] FREES, E.W. Regression Modeling with Actuarial and Financial Applications. Cambridge: Cambridge University Press, 2010. [10] GNEDENKO, B. V. Sur la distribution limite du terme maximum d’ une série aléatoire. Ann. Math. 1943, (44), pp. 423 -–453. [11] GONZÁLES-ESTRADA, E., J.A. VILLASENOR-ALVA. A Goodness-of-Fit Test for Location-Scale Max-Stable Distributions. Communications in Statistics - Simulation and Computation. 2010, 39(3), pp. 557–562. [12] HOSKING, J., J. WALLIS a E. WOOD. Estimation of the Generalized Extreme-Value Distribution by the Method of Probability-Weighted Moments. Technometrics, 1985, 27(3), pp. 251-261. [13] HOSKING, J., J. WALLIS. Parameter and Quantile Estimation for the Generalized Pareto Distribution. Technometrics, 1987, 29(3), pp. 339–349. – 59 – Seznam použité literatury 60 [14] JENKINSON, A. F. The frequency distribution of the annual maximum (or minimum) values of meteorological elements. Q.J.R. Meteorol. Soc. 1955, (81), pp. 158–171. [15] KLUGMAN, S. A., H. H. PANJER a G. E. WILLMOT. Loss models: from data to decisions. 3rd ed. New Jersey: John Wiley & Sons, c2008. ISBN 978-0-470-18781-4. [16] LANDWEHR, J., N. MATALAS a J. WALLIS. Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resour. Res., 1979, (15), pp. 1055 -–1064. [17] MCNEIL, A. J., R. FREY a P. EMBRECHTS. Quantitative Risk Management: Concepts, Techniques and Tools. New Jersey: Princeton University Press, 2005. ISBN 978-0-691-12255-7. [18] PICKANDS J. Statistical Inference Using Extreme Order Statistics. The Annals of Statistics. 1975, 3(1), pp. 119–131. [19] RESNICK, S. I. Extreme Values, Regulas Variation and Point Processes. New York: Springer, 1987. [20] SCARROTT C., A. MACDONALD. A Review of Extreme Value Threshold Estimation and Uncertainty Quantification. REVSTAT - Statistical Journal. 2012, 10(1), pp. 33–60. [21] SMITH, R. L. Maximum likelihood estimation in a class of nonregular cases. Biometrika. 1985, (72), pp. 67–-90. [22] de SOUSA, B. a G. MICHAILIDS. A Diagnostic Plot for Estimating the Tail Index of a Distribution. J. Comput. Graph. Statist. 2004, 13(4), pp. 974–995. [23] VILLASENOR-ALVA, J. A. a E. G. GONZÁLEZ. A bootstrap goodness of fit test for the generalized Pareto distribution. Computational Statistics and Data Analysis. 2009, 53, pp. 3835–3841.