TESTY A ODHADY PARETOVA INDEXU

ROBUST’2004

c JČMF 2004

TESTY A ODHADY PARETOVA INDEXU Jan Picek Klíčová slova: Paret˚ uv index, rozdělení extrémních hodnot, sféra přitažlivosti, Hill˚ uv odhad. Abstrakt: Nechť X1 , X2 , . . . jsou nezávislé stejně rozdělené náhodné veličiny s distribuční funkcí F a nechť Mn = max(X1 , . . . , Xn ). Pro většinu obvyklých distribučních funkcí vhodně standardizovaná maxima Mn konvergují v distribuci k rozdělení extrémních hodnot Gγ . Podle hodnot shape parametru γ rozlišujeme tři základní třídy distribučních funkcí: γ > 0 – Fréchetova třída, γ = 0 Gumbelova a γ < 0 Weibullova. Z hlediska extrémních událostí je především zajímává třída Fréchetova, γ se v tomto kontextu často nazývá Paretovým indexem. V příspěvku se proto budeme zabývat semiparametrickými odhady γ především pro tuto třídu a testy o γ, zvláště se bude jednat o testy hypotézy γ = 0 proti alternativě γ > 0, tj. náhodný výběr je z rozdělení, který patří do Gumbelovy třídy proti alternativě, že rozdělení je z Fréchetovy třídy.

1

Úvod

Nechť X1 , X2 , . . . jsou nezávislé stejně rozdělené náhodné veličiny s distribuční funkcí F . Naše pozornost v tomto článku bude soustředěna na extremální události. Nechť tedy Mn = max(X1 , . . . , Xn ). Zřejmě distribuční funkce Mn je P (Mn ≤ x) = P (X1 ≤ x, . . . , Xn ≤ x) = F n (x) s.j.. Jednoduše je potom možné ukázat, že Mn → xF s.j. pro n → ∞, kde xF := sup{x ∈ IR : F (x) < 1} ≤ ∞. Tato skutečnost nám neposkytne příliš mnoho informace. Pokud se inspirujeme centrální limitní větou, jistě je přirozené se zabývat standardizovanými maximy. Předpokládejme, že m˚ užeme najít posloupnost reálných čísel an > 0 a bn tak, že posloupnost (Mn − bn )/an konverguje v distribuci, t.j. P ((Mn − bn )/an ≤ x) = F n (an x + bn ) → G(x), n → ∞,

(1)

pro nějakou nedegenerovanou d.f. G(x) Jestliže podmínka platí, říkáme, že F je ve sféře přitažlivosti G (domain of attraction) – F ∈ MDA(G). Přirozeně nás patrně napadnou otázky: jak vypadá G, jaké podmínky musí F splňovat, aby F ∈ MDA(G) a jak volit an a bn . Odpověď na tyto základní otázky m˚ užeme najít např. v [2].

276

Jan Picek

Odpověď na první otázku známe už od roku 1928 – Fisherova-Tippettova věta: Jestliže F ∈ MDA(G) potom G je typu jedné z následujících tří distribučních funkcí: 0, x≤0 γ>0 Fréchet Φ1/γ (x) = exp −x−1/γ , x > 0 exp −(−x)1/γ , x ≤ 0 Weibull Ψ1/γ (x) = γ<0 1 x>0 Gumbel Λ(x) = exp (−e−x ) , x ∈ IR. Po vhodné reparametrizaci m˚ užeme tyto tři třídy charakterizovat jediným rozdělením – zobecněnéným rozdělením extrémních hodnot (Generalized Extreme Value Distribution) exp −(1 + γx)−1/γ γ 6= 0 G(x) = Gγ (x) = , exp(−e−x ) γ=0

kde 1 + γx > 0. Hodnota „shapeÿ parametru γ > 0 odpovídá Fréchetově třídě, γ = 0 Gumbelově a γ < 0 Weibullově. Fisherova-Tippettova věta nám pak říká: jestliže vhodně standardizované maxima konvergují v distribuci k nedegenerované limitě, potom limitní rozdělení musí být rozdělení extrémních hodnot. Poznamenejme, že G je určena jednoznačně až na parametr polohy a měřítka. Je možné ukázat, že v podstatě všechny běžně uvažované spojité rozdělení splňují podmínku (1). Než se zaměříme na volbu an a bn připomeňme několik pojm˚ u z klasické teorie extrémních událostí. Funkce h(t) na (0, ∞) je pravidelně se měnící funkce (regularly varying) v ∞ s indexem α ∈ IR (h ∈ Rα ), jestliže lim

x→∞

h(xt) = tα , t > 0. h(x)

Funkce L(t) na (0, ∞) je pomalu se měnící funkce (slowly varying) v ∞ (L ∈ R0 ), jestliže L(xt) lim = 1, t > 0. x→∞ L(x) V oblasti extrémních hodnot se často pracuje s kvantilovou funkcí chvostu 1 −1 U (t) = F 1− = inf{y : F (y) ≥ 1 − 1/t}, t > 0. t Věta 1.1. a) F ∈ MDA(Gγ ) právě když U (tx) − U (t) xγ − 1 = t→∞ a(t) γ lim

pro každé x > 0, a je nějaká kladná funkce a γ ∈ IR, an = a(n), bn = U (n).

277

Testy a odhady Paretova indexu b) F ∈ MDA(Gγ ), γ > 0 právě když lim

t→∞

U (tx) = xγ U (t)

(2)

pro každé x > 0 s γ > 0, tj. U ∈ Rγ (an = U (n)). D˚ ukaz a detaily např. v de Haan L. (1970).

Další a často používá charakterizace Fréchetovy třídy: • F ∈ MDA(Gγ ), γ > 0 právě když 1 − F (x) ∈ R−1/γ , tj. chvost rozdělení F je pravidelně se měnící funkce v ∞ s indexem −1/γ •

1 − F (x) = x−1/γ L(x).

(3)

Statistickou inferenci v extremální statistice m˚ užeme založit na základě limitního rozdělení, tj. na zobecněném rozdělením extrémních hodnot např. pomocí metody maximální věrohodnosti. Ukazuje se, že konvergence je však velmi pomalá, proto je nutné hledat alternativní přístupy. V následujícím textu ukážeme některé možné semiparametrické přístupy.

2

Testy

Případ F ∈ MDA(G0 ) je zajímavý pro mnoho aplikací, které se zabývají extrémy. D˚ uvodem je nejen jednodušší inference založená na Gumbelově sféře přitažlivosti, ale také široká paleta rozdělení s exponenciální chvosty. Jako zástupce jmenujme normální, lognormální a gamma rozdělení. Na druhé straně opravdu extrémní události jsou modelovány pomocí rozdělení z Fréchetovy třídy. Je tedy určitě v praxi užitečné rozhodnout do jaké třídy rozdělení našich dat patří. To znamená uvažovat následující test oboustranné hypotézy (respektive anologický jednostranný test) F ∈ MDA(G0 ) proti alternativě F ∈ MDA(Gγ )γ6=0 .

(4)

Asi nejpoužívanější test pro tuto situaci navrhli Hasofer A.M. and Wang Z. v roce 1992. Najdeme ho implementovaného v řadě softwar˚ u pro statistiku extrémních událostí. Test jako většina semiparametrických postup˚ u je založen na k největších pořádkových statistikách: 2 k k X k − Xn−k+1:n 1X , X := Xn−i+1:n . (5) Wk = k 2 Pk k i=1 (k − 1) i=1 Xn−i+1:n − X k Hasofer a Wang ukázali, že testová statistika Wk má asymptoticky normální rozdělení se střední hodnotou µk a rozptylem σk2 µk =

1 , (k − 1)

σk2 =

4(k − 2) 1 2 (k − 1) (k + 1)(k + 2)

278

Jan Picek

Kritický obor pro oboustrannou alternativu je potom dán následovně |Wk∗ | > u1−α/2 , kde Wk∗ := (Wk − µk )/σk a uε je ε-kvantil normálního rozdělení. Při praktickém provádění testu jistě narazíme na problém, jak zvolit vhodné k. Pokud budeme k zvyšovat, zvýšíme sílu testu, ale na druhé straně zvyšující se podíl k/n má neblahý vliv na chybu I.druhu. Volba se pak stává do jisté míry „alchymiíÿ, nicméně v literatuře existují doporučení, např. Boos navrhuje k/n=0.2 pro 50 √ ≤ n ≤ 500 a k/n = 0.1 pro 500 < n ≤ 5000, Galambos radí volit k = 2 n. Podobný typ testu navrhli C. Neves, J. Picek a M. I. Fraga Alves (2005). Jako testovou statistiku uvažují ∗ Tk,n = 1 k

k P

i=1

Xn:n − Xn−k:n (Xn−i+1:n − Xn−k:n )

− log k .

(6)

∗ Ukázali, že testová statistika Tk,n za nulové hypotézy konverguje k Gumbe−x lovu rozdělení G(x) = exp(−e ) a že test je konzistentní. Nulová hypotéza ∗ je tedy zamítnuta na asymptotické hladině α ∈ (0, 1) jestliže Tk,n < gα/2 ∗ nebo Tk,n > g1−α/2 , kde gε označuje ε-kvantil Gumbelova rozdělení, tj. gε = − log(− log ε).

Jako poslední přístup pro test (4) uveďme poměrně nedávný přístup J. Segerse a J.Teugelse (2001). Vychází z poměru uvažovaném Galtonem (1902): Gn =

Xn:n − Xn−2:n Xn−1:n − Xn−2:n

Náhodný výběr o rozsahu n je rozdělen do m skupin je spočítán poměr (i)

ξi =

Pm

i=1

ni = n. V každé

(i)

Xni :ni − Xni −2:ni (i)

(i)

Xni −1:ni − Xni −2:ni

,

1, · · · , m

Podle Serflinga (1980), Segers a Teugels navrhují užít testovou statistiku Sm

5 = m

m X i=1

!2

T (ξi )

,

T (x) := 1 −

6x , (1 + x)2

a ukazují, že za nulové hypotézy konverguje k χ21 rozdělení pro m → ∞. Nulová hypotéza je tedy zamítnuta na asymptotické hladině α, je-li Sm > χ21 (1 − α), kde χ21 (ε) označuje ε-kvantil χ2 rozdělení s 1 st. vol.

(7)

279

Testy a odhady Paretova indexu

2.1

Numerická ilustrace

0.0

0.2

0.4

0.6

0.8

1.0

Zkusme ilustrovat chování výše uvedených test˚ u na simulovaných datech a na jednom reálném příkladu. Nejprve jsme uvažovali platnost nulové hypotézy (4), tj. jako zástupce z Gumbelovy sféry přitažlivosti jsme zvolili Gumbelovo rozdělení F (x) = exp(−e−x ). Z tohoto rozdělení jsme vygenerovali 1000 × výběr o rozsahu 1000 a provedli výše uvedené testy. Na obr. 1 jsou zobrazeny výsledky ve formě relativního počtu zamítnutí nulové hypotézy na hladině α = 0.05. Testy (5) a (6) byly provedeny pro k = 2, . . . , 999 (počet použitých nejvyšších pořádkových statistik). Test (7) byl konstruován tak, že výběr byl rozdělen do 50 (=m) blok˚ u o rozsahu 20. Obr. 1 vlastně ilustruje odhad chyby prvního druhu. Je vidět, že odhad této chyby pro test (7) je prakticky 0.05, pokud přijmeme výše zmiňovaná doporučení pro volbu k, potom testy (5) a (6) mají odhad také blízký 0.05. Nicméně se zdá, že test (6) dovolí volit větší rozsah k aniž by to mělo výrazný vliv na chybu I. druhu. Testovali jsme i jiná rozdělení z Gumbelovy sféry přitažlivosti i pro jiné rozsahy, charakter křivek byl podobný s jedinou výjimkou a to exponenciálním rozdělením, pro které odhad chyby prvního druhu byl stabilní (blízko hodnoty 0.05) prakticky pro všechna možná k.

0

200

400

600

800

1000

k

Obrázek 1: Relativní počet zamítnutí H0 na hladině α = 0.05 pro Gumbelovo ∗ rozdělení, Tk,n (plná čára), Wk (tečkovaně), S50 (čerchovaně). Jako další zástupce pro ilustraci bylo zvoleno zobecněné Paretovo rozdělení 1 Fγ (x) := 1 + log Gγ (x) = 1 − (1 + γx)− γ x≥0 jestliže γ ≥ 0 pro 0 ≤ x ≤ − γ1 jestliže γ < 0

280

Jan Picek

Toto rozdělení závisí na parametru γ. Podle jeho hodnoty patří rozdělení do jedné z uvažovaných tříd. Opět byl 1000 krát generován výběr o rozsahu 1000 pro hodnoty γ = −2.0, -1.5, -1.0, -0.5, -0.25, -0.1, -0.01, 0.01, 0.1, 0.25, 0.5, 1.0, 1.5, 2.0. Pokud se opět zajímáme o relativní počet zamítnutí nulové hypotézy, pak v tomto kontextu dostáváme představu o síle test˚ u. Na obr. 2 vidíme srovnání pro všechny tři testy v závislosti na γ data. Testy (5) a (6) byly provedeny pro k = 150, test (7) s m = 50. Ve všech třech případech tak bylo použito 150 hodnot (i když ne nutně stejných).

1

power

0.8

0.6

0.4

0.2

0 -2

-1

0 gamma

1

2

∗ Obrázek 2: Síla testu: T150,n (plná), W150 (čerchovaná), S50 (tečkovaná) na hladině α = 0.05 pro zobecněné Pareto (γ = −2.0, -1.5, -1.0, -0.5, -0.25, -0.1, -0.01, 0.01, 0.1, 0.25, 0.5, 1.0, 1.5, 2.0), rozsah n = 1000.

Vidíme, že z hlediska síly testu se nejlépe chová test (5), trochu h˚ uře (6) a nejslabší je test (7). Ten byl nejslabší ve všech případech, které jsme zkoumali. Testy (5) a (6) se příliš nelišily a záviselo na konkrétní volbě k, rozdělení a rozsahu. Dokladem toho m˚ uže být např. obr. 3, který zobrazuje závislost „síly testuÿ na volbě k pro zobecněné Pareto rozdělení s γ = 1.0 Co se týče asymptotických vlastností a předpoklad˚ u všechny tři testy jsou rovnocenné, na druhou stranu vidíme, že pokud máme i poměrně velký rozsah dat, rozdíly najít m˚ užeme. Nejslabším testem se zdá být do jisté míry (7). Test (5) je v praxi patrně nejpoužívanější, ale zdá se, že (6) je plně srovnatelný. Podívejme se též na testy na reálných datech. V poslední době se vede diskuse, že počasí nabývá extrémního chování. Jedním z mnoha charakteristik tohoto chování počasí mohou být např. extrémní srážky. V České Republice jsou k dispozici data na řadě stanic od roku 1961. Extrémní srážky m˚ užeme

281


1

power

0.8

0.6

0.4

0.2

0 0

50

100 k

150

200

50

100

150

200

250

300

∗ Obrázek 3: Síla testu: Tk,n (plná čára), Wk (čerchovaná), S20 (tečkovaná) na hladině α = 0.05 pro v závislosti na k, rozsah n = 200 pro zobecněné Pareto rozdělení s γ = 1.0 (vpravo).

1961

1970

1980

1990

2000

Obrázek 4: Maximální třídenní úhrny srážek v letech 1961-2000 ve Valašském Meziříčí. třeba charakterizovat maximálními třídenními úhrny srážek v daném roce (takovéto data měl autor k dispozici). Na obr. 4 vidíme tyto data pro stanici ve Valašském Meziříčí. Velmi dobře je vidět výjimečný rok 1997, který přinesl velké záplavy na Moravě. Je otázkou pro další statistické úvahy, jaký základní model je pro tuto veličinu

282

Jan Picek

-4

-2

0

2

4

6

(maximálními třídenními úhrny srážek v daném roce) vhodný, tj. Gumbelova nebo Fréchetova třída. Výsledky test˚ u jsou graficky zobrazeny na obr. 5, kde vodorovné čáry odpovídají příslušným 97.5%-ním kvantil˚ um pro oboustranný test. Vidíme, že zamítnutí nulové hypotézy je velmi problematické, zamítáme pouze pro větší hodnoty k a to hlavně testem (6), viděli jsme ze simulací, že větší k nedávají dobré výsledky co se týče platnosti nulové hypotézy. Hlavním problémem tu je však velmi malý počet pozorování (n = 40), který je v aplikacích týkajících se extrému nedostatečný, ale bohužel v praxi častý.

0

10

20

30

40

∗ Obrázek 5: Srážky ve Valašském Meziříčí: Hodnoty Tk,40 (plná), Wk (tečkovaná), S8 (čerchovaná). Vodorovné linky označují příslušné kvantily odpovídající α = 0.05.

Další testy, které byly v poslední době konstruovány, uvažují hypotézy o hodnotách parametru γ (jak „těžkéÿ jsou těžké chvosty) pro F ∈ MDA(Gγ ), γ > 0, viz [11], [16]. O těchto testech bylo referováno na Robustu 2000. Pokud se budeme zabývat úvahami o hodnotách γ, pak mnohem bohatší je literatura věnovaná odhad˚ um. Proto následující část tohoto příspěvku věnujeme právě jim.

3

Odhady

Připomeňme, že vycházíme z náhodného výběru X1 , X2 , . . . z rozdělení s neznámou distribuční funkcí F . Pokud F ∈ MDA(Gγ ), γ > 0, pak patrně nejznámější odhadem γ je Hill˚ uv odhad z roku 1975 [8]: Hn (k) =

k−1 1X log X(n−i:n) − log X(n−k:n) . k i=0

(8)

283


Ukažme návrh jedné z možných cest jeho odvození. Vyjděme z charakterizace Fréchetovy třídy (2): lim

t→∞

U (tx) = xγ , kde U (t)

U (t) = F −1 (1 − 1/t).

Po zlogaritmování dostaneme limt→∞ log U (t/x) − log U (t) = −γ log x. Výběrová verze kvantilové funkce chvostu U je Un (1/x) = Fn−1 (1 − x) = n Xn(1−x),n , tj. Un ( nk ) = Xn−k,n a Un ( kx ) = Xn−kx,n . Tedy pro 0 < x < 1 je log Xn−kx,n − log Xn−k,n = −γ log x. Poté integrujme γ = −γ

Z

1

log x dx = lim

t→∞

0

Z

0

1

{log U (t/x) − log U (t)} dx.

Dostaneme tak možný odhad γ Hn (k) =

Z

0

=

1

(log Xn−kx,n − log Xn−k,n ) dx

k−1 1X log X(n−i:n) − log X(n−k:n) k i=0

Hill˚ uv odhad je konzistentní, tvrzení najdeme např. v [13]. Věta 3.1. Je-li F ∈ MDA(Gγ ), γ > 0, potom Hn (k) → γ v pravděpodobnosti, k = k(n) → ∞, k(n)/n → 0(n → ∞). Pokud nás zajímá asymptotické rozdělení odhadu, musíme klást další podmínky na distribuční funkci, abychom byli schopní ho odvodit. Nejčastěji se uvažuje následující podmínka (regular variation of second order): Nechť existuje A(t) funkce konstantního znaménka a parametr ρ lim

t→∞

U(tx) U(t)

− xγ

A(t)

= xγ ·

xρ − 1 ρ

(9)

pro všechna x > 0. Věta 3.2. Nechť√podmínka (9) platí a nechť posloupnost k = k(n) je taková, že k(n) → ∞ a kA(n/k) → 0, potom √ k(Hn (k) − γ)

je asymptoticky normální s nulovou střední hodnotou a rozptylem γ 2 .

284

Jan Picek

Pokud uvažujeme F ∈ MDA(Gγ ), γ libovolné, pak lze „analogicky odvoditÿ Momentový odhad [1] (1)

M (k) = 1 + M (k) kde M (k)(j) =

1 + 2

!−1 2 M (k)(1) −1 , M (k)(2)

(10)

k j 1X log X(N n−i+1:N n) − log X(N n−k:N n) . k i=1

Alternativou momentového odhadu je Pickands˚ uv odhad [17] 1 XN n−k+1:N n − XN n−2k+1:N n P (k) = log . log 2 XN n−2k+1:N n − XN n−4k+1:N n

(11)

Výše uvedené odhady jsou patrně nejznámější, v literatuře existuje obrovské množství dalších odhad˚ u: r˚ uzná zobecnění Hillova odhadu, odhady založené na parametru druhého řádu ρ, viz (9) a mnoho a mnoho dalších alternativ. Uveďme alespoň jeden příklad, který vychází z (9) a uvažuje, že ρ = −1. Navrhli ho Gomes a Martin v roce 2002, viz [9]. k 1X Ui − GM (k) = k i=1

k 1X iUi k i=1

! P k

i=1 (2i

Pk

i=1

− k − 1)Ui

i(2i − k − 1)Ui

XN n−i+1:N n Ui = i log , XN n−i:N n

,

(12)

Stejně jako u test˚ u je problém volby k, lze řešit podobnými doporučeními nebo se uvažují postupy založené na bootstrapu - viz např. [10]. Pokud se podíváme do domácích luh˚ u a háj˚ u, tak i tady najdeme příspěvek ke konstrukci odhad˚ u parametru γ za podmínky F ∈ MDA(Gγ ), γ > 0. Tyto odhady nejsou založeny přímo na pořádkových statistikách na rozdíl od předcházejících. Vychází se opět z určité charakterizace Fréchetovy třídy: − log(1 − F (a)) = 1. a→∞ m log a lim

(13)

Aplikací l’Hospitalova pravidla z von Mises podmínek (viz Embrechts a kol., Kap. 3), dostaneme, že 1 − F (x) = x−m L(x), což je charakterizace uvedená v (3). Platí i opačná implikace. Principem spočívá v rozdělení výběru do skupin, v každé je spočtena nějaká jednoduchá statistika. Výsledný odhad je konstruován na základě empirické distribuční funkce sledované statistiky.

285


O prvním typu odhadu referovala na Robustu 2000 A. Fialová: Rozdělíme pozorování do N nepřekrývajících se výběr˚ u rozsahu n a určíme zde pr˚ uměry (1) (N ) ¯ ¯ (Xn , . . . , Xn ). Dostaneme pak náhodný výběr z rozdělení s distribuční ) 1 PN ¯ n ≤ x). Označíme Fb (N ¯ (j) ≤ x] funkcí FX¯n (x) = IP (X ¯ n (x) = N j=1 I[Xn X ¯ n(1) , . . . , X ¯ n(N ) ). Vyberme poempirickou distribuční funkcí založenou na (X

sloupnost {aN }∞ N =1 , aN → ∞ pro N → ∞ ve tvaru aN = N ným δ ∈ (0, 1).

1−δ m0

, s pev-

Odhad parametru m = 1/γ je potom (N )

(N )

m bN = m ˜ N (aN )I[0 < FbX¯ n (aN ) < 1] + m0 I[FbX¯ n (aN ) = 0 nebo 1],

kde

m ˜ N (a) =

(N ) − log 1 − FbX¯n (a) log a

,

(14)

a > 0.

Odhad (14) je konzistentní a jeho asymptotické rozdělení je normální, viz následující věty. Věta 3.3. Nechť {X1 , X2 , . . .} je posloupnost nezávislých stejně rozdělených náhodných veličin s distribuční funkcí F ∈ MDA(Gγ ), γ > 0 a hustotou f (x) = 0 pro x < 0 a 0 < f (x) < ∞ for x ≥ Kf ≥ 0. Nechť m b N je odhad m. Potom m b N → m s pravděpodobností 1, pro N → ∞. Věta 3.4. Za podmínek předcházející věty posloupnost 1 1 − FX¯ n (aN ) 2 1 log L∗ (aN ) N 2 log aN m bN − m + FX¯n (aN ) log aN

je asymptoticky normální N → ∞.

D˚ ukazy obou vět lze nalézt v [4]. Na rozdíl od Hillova (a i dalších výše zmíněných) odhad˚ u je asymptotické rozdělení odvozeno za mnohem slabších předpoklad˚ u. Bohužel výsledek věty 3.4 obsahuje pomalu měnící se funkci L∗ , kterou zpravidla neznáme, není proto možné jednoduše výsledku využít např. pro konstrukci intervalových odhad˚ u. Pro tento odhad musíme zvolit δ, což je vlastně podobná úloha jako je určení vhodného k pro předcházející odhady, navíc je však nutné zvolit m0 , což vyžaduje nějakou počáteční informaci o tom, jak chvost rozdělení m˚ uže být „těžkýÿ. To poněkud omezuje užití odhadu pro praktické problémy. Ilustrujme na simulovaných datech na chování odhadu právě v závislosti na volbě δ a m0 . Jako model dat použijeme Paretovo rozdělení, které je jedním z typicky používaných rozdělení pro popis extrémních událostí: 1/γ 1 F (x) = 1 − , x≥0 (15) 1+x

286

Jan Picek

1.2 1.0 0.8 0.6

estimate

0.4

0.6 0.4

estimate

0.8

1.0

1.2

Konkrétně byla simulace provedena pro Paretovo rozdělení s γ = 1, což je i hodnota, kterou chceme odhadnout. Výsledek m˚ užeme vidět na obr. 6, z kterého vyplývá, že pokud je zhruba δ < 0.4 je odhad poměrně stabilní a rozumný. Pro velké hodnoty δ odhad naprosto selhává. Zároveň je vidět, že čím horší máme apriorní informaci o správné hodnotě γ = 1/m, tím dostaneme horší výsledek.

0.0

0.2

99% 75% 50% 25% 1%

0.0

0.0

0.2

99% 75% 50% 25% 1% 0.2

0.4

0.6 delta

0.8

1.0

1

2

3

4

5

6

m_0

Obrázek 6: Závislost odhadu v 1000 simulovaných výběrech Paretova rozdělení s γ = 1 na parametru δ pro dané m0 = 1.5 (vlevo) závislost hodnot odhadu na parametru m0 pro δ = 0.1 (vpravo). Uvedeny jsou medián 1, 25, 75 a 99 percentily. Jurečková, Picek (2004) navrhli odhad vycházející z postup˚ u pro testování hypotézy o hodnotách parametru γ pro F ∈ MDA(Gγ ), γ > 0. Krátká poznámka o nich byla v předcházející kapitole. Invertováním těchto test˚ u (v duchu zp˚ usobu, který užil Hodges a Lehmann v roce 1963) dostaneme odhad 1 + − MN = (MN + MN ), 2 kde − MN = sup{s : 1 − FˆN∗ (aN,s )) < N −(1−δ) }, + MN = inf{s : 1 − FˆN∗ (aN,s )) > N −(1−δ) }. (1)

(N )

X(n) , . . . , X(n) jsou odpovídající výběrová maxima N skupin o rozsahu n, které vznikly rozdělením p˚ uvodního náhodného výběru. Jako FˆN∗ označujeme empirickou distribuční funkci odpovídající výběrovým maxim˚ um, aN,m = (nN 1−δ )1/m , kde 0 < δ < 21 je konstanta. Ilustrujme podobně jako u předcházejícího odhadu chování v závislosti na volbě δ. Jako model dat tentokráte použijeme Burrovo rozdělení, které je dalším typicky používaným rozdělením pro popis extrémních událostí: α 1 F (x) = 1 − , x≥0 (16) 1 + x1/γ

Konkrétně byla simulace provedena pro Burrovo rozdělení s γ = 1, α = 1, jednička je opět hodnota, kterou chceme odhadnout. Výsledek m˚ užeme vidět

287


1.3

na obr. 7, z kterého vyplývá, že lepší výsledek dostaneme, pokud je δ blízké 0.5. Neplatí to obecně, pro jiná rozdělení to m˚ uže dopadne úplně opačně. Na druhou stranu volba δ není tak problematická jako volba k u Hillova odhadu, viz obr. 8, kde na stejných datech je spočítán Hill˚ uv odhad.

1.1 0.9

1.0

odhad

1.2

5% 25% 50% 75% 95%

0.0

0.1

0.2

0.3

0.4

0.5

delta

Obrázek 7: Závislost odhadu v 1000 simulovaných výběrech Burrova rozdělení s γ = 1, α = 1 na parametru δ. Uvedeny jsou medián 5, 25, 75 a 95 percentily.

1.0 0.4

0.6

0.8

odhad

1.2

1.4

1.6

5% 25% 50% 75% 95%

0

100

200

300

400

500

k

Obrázek 8: Hill˚ uv odhad v 1000 simulovaných výběrech o rozsahu 1000 v závislosti na k pro Burrovo rozdělení s γ = 1, α = 1. Uvedeny jsou medián 5, 25, 75 a 95 percentily.

288

Jan Picek

Jurečková a Picek ukázali v [12], že odhad je silně konzistentní. Asymptotickou normalitu odvodil Omelka [15]. Odhad (16) potřebuje pouze volbu δ, což ho činí použitelnějším než odhad (14). I simulace dávájí poměrně dobré výsledky - viz dále, přesto však musíme být v praktických aplikacích velmi opatrní. Oba odhady nejsou invariatní vzhledem ke změně měřítka na rozdíl od odhadu Hillova (8), Pickandsonova (11), momentového (10) i (12). Všechny zmíněné odhady nejsou invaritní vzhledem ke změně polohy. Při mechanickém použití odhad˚ u to potom m˚ uže vést k „zajímavýmÿ výsledk˚ um. Byly proto uvažovány některé modifikace Hillova odhadu, viz např. [3]. Nějaké poznámky, jak se s naznačeným problémem vypořádat pro odhad (16) učinil Omelka [15].

3.1

Numerická ilustrace

V této části zkusíme porovnat zmiňované odhady na simulovaných datech. Jako výchozí model použijeme dříve zmiňovaná rozdělení: Paretovo, Burrovo a zobecněné Paretovo. U všech tří rozdělení zvolíme shape parametr (γ = 1/m = 1) tak, aby chvosty byly „stejně těžkéÿ a mohli tak sledovat vliv rozdělení. U zobecněného Paretova zvolíme ještě další dvě hodnoty γ: 1/3 – lehčí a 2 – těžší chvost. Z daného rozdělení jsme vygenerovali 1000 × výběr o rozsahu 1000 a provedli výše zmíněné odhady. Odhady (8), (11), (10) a (12) jsme spočítali pro k = 2, . . . , 998. Pro odhady (14) a (16) jsme provedli rozdělení do 200 skupin (= N ) po 5 hodnotách (= n) a spočítali odhad pro δ = 0.01, . . . , 0.50 s krokem 0.01, navíc pro (14) za m0 jsme zvolili „skutečnéÿ 1/γ + 1. Za účelem porovnání jsme pro každé k, respektive δ spočetli střední kvadratickou chybu (MSE) a vybrali takovou hodnotu k (δ), kdy je MSE minimální a spočítali nějaké výběrové charakteristiky z tisíce získaných hodnot odhad˚ u. Výsledky najdeme v tabulce 1. Odhad (14) je v ní označen jako FJP a (16) jako JP. Tučně je zvýrazněna pro dané rozdělení minimální MSE mezi odhady. M˚ užeme si všimnout, že pro opravdu těžké chvosty, tj. pro všechny případy kromě zobecněného Paretova rozdělení s γ = 1/3, dávají všechny odhady v pr˚ uměru rozumné výsledky. Nejslabší se přesto zdá být odhad (14) a protože byly už některé výhrady diskutovány dříve, nelze ho doporučit pro praktické úlohy. Naopak odhad (16) je srovnatelný s klasickými, navíc pro lehčí chvosty dává často rozumnější výsledky než klasické odhady. Zdá se tedy, že s ním lze pracovat minimálně jako vhodnou alternativou. Z tabulky je dále vidět a další simulace pro jiné případy a rozdělení to jen potvrzují, že index lehčích chvost˚ u se odhaduje mnohem h˚ uře. Překvapující je výsledek Pickandsova odhadu (alespoň pro autora tohoto příspěvku), protože tento odhad by měl fungovat pro odhad nejen ve Fréchetově třídě, ale i pro Gumbelovu a Weibullovu sféru přitažlivosti, tedy i pro „lehkéÿ chvosty.

289


rozdělení

metoda

k, δ

MSE

pr˚ uměr

medián

rozptyl

Pareto γ=1

Hill Moment Pickands Gomes FJP JP

k = 998 k = 998 k = 985 k = 997 δ = 0.15 δ = 0.49

0.0010 0.0023 0.0221 0.0044 0.0123 0.0147

1.0003 1.0053 1.0177 1.0016 0.9542 1.0435

0.9984 1.0033 0.9967 0.9968 0.9371 1.0427

0.0010 0.0022 0.0218 0.0044 0.0102 0.0128

Burr α=1 γ=1


k = 112 k = 257 k = 985 k = 998 δ = 0.22 δ = 0.49

0.0098 0.0101 0.0221 0.0012 0.0111 0.0047

0.9517 0.9478 1.0177 1.0007 0.9574 1.0226

0.9489 0.9383 0.9967 0.9989 0.9402 1.0192

0.0075 0.0074 0.0218 0.0012 0.0093 0.0042

zobec. Pareto γ=2 β=1


k = 310 k = 367 k = 993 k = 482 δ = 0.01 δ = 0.45

0.0010 0.0010 0.0020 0.0025 0.0084 0.0042

0.4847 0.4880 0.5030 0.5227 0.4177 0.5519

0.4841 0.4863 0.4997 0.5210 0.4123 0.5491

0.0007 0.0009 0.0020 0.0020 0.0016 0.0015

zobec. Pareto γ=1 β=1


k = 112 k = 257 k = 985 k = 998 δ = 0.22 δ = 0.49

0.0098 0.0101 0.0221 0.0012 0.0111 0.0047

0.9517 0.9478 1.0177 1.0007 0.9574 1.0226

0.9489 0.9383 0.9967 0.9989 0.9402 1.0192

0.0075 0.0074 0.0218 0.0012 0.0093 0.0042

zobec. Pareto γ = 1/3 β=1


k = 23 k = 257 k = 890 k = 102 δ = 0.01 δ = 0.11

0.5527 0.5037 16.1112 0.4795 0.2869 0.1166

2.4329 2.5140 3.6237 2.4276 2.5618 2.8500

2.3598 2.4248 3.0364 2.4020 2.5565 2.8384

0.2314 0.2678 15.7379 0.1520 0.0949 0.0942

Tabulka 1: Výběrové charakteristiky odhad˚ u Paretova indexu pro minimální MSE při 1000 opakování generování dat rozsahu 1000 pro r˚ uzná rozdělení.

Reference [1] Dekkers A.L.M, Einmahl J.H.J., de Haan L. (1989). A moment estimator for the index of an extreme value distribution. Ann. Statist. 17, 1833-1855. [2] Embrechts P., Klüppelberg C., Mikosch T. (1997). Modelling extremal events for insurance and finance. Springer-Verlag, Berlin.

290

Jan Picek

[3] Fraga Alves M.I. (2001). A location invariant Hill-type estimator. Extremes, 4 (2), 165 – 183. [4] Fialová A., Jurečková J., Picek J. (2004). Estimation of tail index based on sample mean. Revstat, 2, 75 – 99. [5] de Haan L. (1970). On regular variation and its application to the weak convergence of sample extremes. Mathematical Centre Tract 32, Amsterdam. [6] de Hann L., Stadtmüller U. (1996). Generalized regular variation of second order. J.Austral.Math.Soc. (A) 61, 381 – 395. [7] Hasofer A.M., Wang Z. (1992). A test for extreme value domain of attraction. JASA, 87, 171 – 177. [8] Hill B.M. (1975). A simple general approach to inference about the tail of a distribution. Ann. Statist. 3, 1163 – 1174. [9] Gomes M.I., Martins M.J. (2002). Asymptotically unbiased estimators of the tail index based on the external estimation of the second order parameter. Extremes 5 (1), 5 – 31. [10] Gomez I., Oliviera O. (2001). The bootstrap methodology in statistics of extremes-choice of the optimal sample fraction. Extremes 4 (4), 331 – 358. [11] Jurečková J., Picek J. (2001). A class of tests on the tail index. Extremes, 4,(2), 165 – 183. [12] Jurečková J., Picek J. (2004). Estimates of the tail index based on nonparametric tests. Theory and Applications of Recent Robust Methods, Birkhauser, Basel, 141 – 152. [13] Mason D.M. (1982). Laws of large numbers for sums of extreme values. Ann. Probab. 10, 754 – 764. [14] Neves C., Picek J., Fraga Alves M. I. (2005). The contribution of the maximum to the sum of excesses for testing max-domains of attractions. J. Statist. Planning Infer., v tisku. [15] Omelka M. (2005). Asymptotic normality of the estimates of the tail index based on nonparametric tests. Zasláno. [16] Picek J., Jurečková J. (2001). A class of tests on the tail index using the modified extreme regression quantiles. Sborník konference ROBUST’00 (J.Antoch, G.Dohnal, eds.), JČMF Praha, 217 – 226. [17] Pickands J. (1975). Statistical inference using extreme order statistics. Ann.Statist. 3, 119 – 131. [18] Segers J., Teugels J. (2001). Testing the Gumbel hypothesis by Galton’s ratio. Extremes, 3:3, 291 – 303. Poděkování: Příspěvek vznikl za podpory Grantové agentury AV ČR – projekt B3042303 a výzkumného záměru MSM4674788501. Adresa: J. Picek, Katedra aplikované matematiky, Technická univerzita v Liberci, Hálkova 6, 461 17 Liberec E-mail : [email protected]

TESTY A ODHADY PARETOVA INDEXU

Recommend Documents