Inhoudsopgave
Inhoudsopgave blz. 1
Hoofdstuk 1: Probleemomschrijving blz. 2
Hoofdstuk 2: Extreme waarden theorie blz. 4
Hoofdstuk 3: Schatters voor de GEV-verdeling en het POT-model blz. 9
Hoofdstuk 4: Copula’s blz. 17
Hoofdstuk 5: Modelleren met Copula’s blz. 20
Hoofdstuk 6: Data en Resultaten blz. 26
Hoofdstuk 7: Conclusies blz. 44
Bibliografie blz. 47
Hoofdstuk 1: Probleemomschrijving
Deze scriptie gaat over extremen in het weer. De laatste maanden is er veel aandacht voor het klimaat. Veel mensen maken zich zorgen over temperatuurstijging en andere gevolgen van klimaatsverandering. Mede ge¨ınspireerd door deze aandacht ben ik begonnen aan deze scriptie. Het belangrijkste doel van deze scriptie is echter niet om te onderzoeken of er daadwerkelijk een stijging van de temperatuur heeft plaatsgevonden of dat er meer buien zijn dan vroeger, hoewel dit wel terloops ter sprake zal komen. De focus in dit onderzoek zal liggen in het modelleren van weersextremen en vervolgens afzonderlijke elementen uit het weer met elkaar in verband te brengen. Hiermee probeer ik een antwoord te kunnen geven op vragen als: Hoe groot is de kans dat er in Nederland op een willekeurige dag 100 mm water valt? Is het waar dat als het vandaag regent, de kans groot is dat de wind dan ook sterker is dan wanneer het niet regent? Ik zal aan het eind van de scriptie proberen antwoord te geven op deze vragen.
De techniek die ik in deze scriptie gebruik kan ook worden toegepast op andere belangrijke onderwerpen. Zo is het voor een verzekeraar interessant te weten hoe (extreme) schade-claims in zijn portefeuille gemodelleerd kunnen worden. Hoe groot is de kans dat de totale schade in een jaar meer dan tien miljoen bedraagt? Als deze verzekeraar meerdere por-tefeuilles heeft en meerdere polissen aanbiedt, is het voor hem belangrijk om te weten of en hoe de schadeclaims verband met elkaar houden. Hoe groot is de kans dat zowel in de ene portefeuille als in de andere een grote schadeclaim plaatsvindt? Of omgekeerd: Welk bedrag aan schade van alle portefeuilles samen wordt met 95% kans niet overschreden? Deze bekende Value-at-Risk benadering kan in potentie met behulp van technieken uit deze scriptie bepaald worden.
Het onderzoek naar de extremen en verbanden in het weer zal ik doen aan de hand van een onderzoeksvraag en vier deelvragen. De onderzoeksvraag, die ik probeer te beantwoorden, luidt:
Hoe kunnen de extremen op het gebied van zonneschijn, temperatuur, luchtdruk, windsnel-heid en neerslag gemodelleerd worden en wat is de samenhang tussen deze elementen uit het weer?
Deze vraag zal ik beantwoorden aan de hand van de volgende deelvragen:
1. Hoe worden extremen in het algemeen gemodelleerd?
2. Kan de extreme waarden theorie op de weerselementen worden toegepast?
3. Hoe kunnen afhankelijke reeksen gemodelleerd worden?
Het antwoord op de eerste deelvraag zal grotendeels gegeven kunnen worden na de be-spreking van de theorie in hoofdstuk twee over extreme waarden theorie. Twee belangrijke onderdelen van deze theorie vormen het Peak Over Threshold (POT) model en de Genera-lised Extreme Value (GEV) verdeling. Het Peak Over Threshold model veronderstelt dat overschrijdingen boven een bepaalde grens verdeeld zijn volgens de gegeneraliseerde Pareto verdeling. De keuze van deze grens is een moeilijke exercitie. In hoofdstuk 3 zal ik ingaan op de manier waarop goede POT-modellen gemaakt kunnen worden en hoe gebruik gemaakt kan worden van de GEV-verdeling. Ik introduceer daar verschillende schatters voor beide modellen en bespreek de verschillende eigenschappen. Aan de hand van dit hoofdstuk kan ik beginnen met het modelleren van de verschillende elementen uit het weer.
Hoofdstuk 2: Extreme waarden theorie
In dit hoofdstuk geef ik een overzicht van de belangrijkste en voor mijn onderzoek meest relevante resultaten uit de Extreme Waarden Theorie. Ik introduceer twee belangrijke verdelingen, de Generalised Extreme Value Distribution (GEV) en de Generalised Pareto Distribution (GPD) verdeling. Dit zijn twee veel gebruikte verdelingen om extremen te modelleren. Daarnaast bespreek ik het Peak Over Threshold (POT) model, een bekend model om extremen te modelleren. Dit model modelleert het gedrag van waarnemingen die boven (onder) een bepaalde grens liggen.
De basis van de klassieke extreme waarden theorie is een stelling over de limietverdeling van maxima, de stelling van Fisher-Tipett, zie [1]:
Laat X = (X1, X2, . . . , Xn)n∈N een rij van indentiek onafhankelijk verdeelde stochastische
variabelen. Definieer Mn := max{X1, X2, . . . , Xn}. Stel er zijn constantes cn > 0, dn ∈ R
en er is een niet gedegenereerde verdeling H zodat geldt: c−1n (Mn− dn)
d
→ H, dan behoort H tot een van de volgende typen verdelingen:
Fr¨echet : Φα(x) = 0, x ≤ 0 exp(−x−α), x > 0 α > 0. Weibull : Ωα(x) = exp(−(−x)α), x ≤ 0 1, x > 0 α > 0.
Gumbel : Λ(x) = exp(− exp−x), x ∈ R.
Deze stelling beweert dat ´als het genormaliseerde maximum van een reeks identiek en
onafhankelijk verdeelde stochasten convergeert in verdeling, dan convergeert het
genor-maliseerde maximum naar een Fr´echet-, Weibull- of Gumbelverdeling. Verdelingen van
het type Fr´echet hebben zogenaamde ’dikke staarten’. Voorbeelden van dit type zijn de
In de extreme waarden theorie is veel geschreven over de Generalised Extreme Value (GEV) verdeling, zie ondermeer [2]. Deze verdeling is de generalisatie van de drie eer-der ge¨ıntroduceerde extreme waarden verdelingen en is als volgt gedefinieerd:
H(x) :=
exp(−(1 + ξx)−1/ξ) als ξ 6= 0,
exp(− exp(−x)) als ξ = 0,
waarbij geldt : 1 + ξx > 0.
De relatie tusen de GEV en de drie verdelingen komt tot stand door middel van de para-meter ξ. De volgende relatie geldt:
Hξ = Φ1/ξ als ξ > 0, Λ als ξ = 0, Ω−1/ξ als ξ < 0.
De relatie tussen parameter α in de Fr´echet of Weibull verdeling en parameter ξ van de
GEV-verdeling volgt direct uit de bovenstaande notatie. Een ander begrip uit de extreme waarden theorie is het maximum attractiedomein. Dit definieer ik als volgt:
Laat X = (X1, X2, . . . , Xn)n∈N onafhankelijk en identiek verdeeld. Definieer nu
Mn := max{X1, X2, . . . , Xn} en stel er bestaan constantes cn > 0 en dn ∈ R zodanig dat
geldt:
c−1n (Mn− dn) d
→ H. (1)
In dat geval noteer ik: X ∈ M DA(H).
Met behulp van het begrip maximum attractiedomein kan het Peak-Over-Threshold model worden ge¨ıntroduceerd. Definieer allereerst de Generalised Pareto Distribution (GPD) door:
G(x) := 1 − (1 + ξ)
−1/ξx als ξ 6= 0,
1 − exp(−x) als ξ = 0.
De basis van het POT-model is nu de volgende benadering uit Pickands [3]:
Laat X ∈ MDA(H), waarbij H is de GEV-verdeling. De volgende benadering geldt nu: P (X < x + u|X > u) ≈ G(x),
voor een voldoende grote u en voor x > 0. Naarmate u kleiner is, wordt de benadering minder precies. Omgekeerd, naarmate u groter wordt, is de benadering preciezer.
Deze benadering betekent dat zolang u groot genoeg is, dat overschrijdingen over u bij be-nadering een GPD-verdeling volgen. Op grond van data kunnen vervolgens de parameters van de GPD geschat worden. Hoe groot u precies moet zijn om een toepasbaar resultaat te krijgen kan alleen op grond van de praktijk bepaald worden, dit betekent een lastige keuze. De keus is een afweging tussen de kwaliteit van de benadering en de kwaliteit van de schat-ters voor de GPD-parameschat-ters. Immers hoe hoger de grens u, hoe minder overschrijdingen en daarom minder data op grond waarvan de parameters geschat kunnen worden. Daaren-tegen, hoe lager u hoe minder goed de benadering van de verdeling van de overschrijdingen door de GPD zijn. In de literatuur zijn verschillende manieren bekend om een keus voor u te maken, maar er is hiervoor geen procedure die altijd werkt en als beste beschouwd kan worden. Ik geef hier een vrij eenvoudige procedure die vaak gebruikt wordt en in eerdere artikelen goede resultaten heeft opgeleverd. Andere methodes zijn te ingewikkeld of tijd-rovend en niet per se beter. Ik introduceer eerst de gemiddelde overschrijdingsfunctie: Laat X een willekeurige stochast en laat u ∈ R. Definieer de gemiddelde overschrijdings-functie e(u) door:
e(u) := E[X − u|X > u]. Het volgende kan nu worden aangetoond:
Laat X = (X1, X2, . . . , Xn)n∈N onafhankelijk en identiek GPD-verdeeld, onder deze
voor-waarde geldt:
e(u) = β + ξu
1 − ξ ,
voor u in het domein van de GPD-verdeling en voor ξ < 1 en waarbij β en ξ zijn de parameters van de GPD verdeling.
Ik introduceer nu de empirische gemiddelde overschrijdingsfunctie eemp(u), dat is het ge-middelde van alle waarnemingen die boven deze grens komen vermindert met u oftewel:
eemp(u) = P
i:Xi>u
Xi−u
Nu , waarbij Nu = card{i : Xi > u}.
De grens u moet nu zo gekozen worden dat de empirische gemiddelde overschrijdingsfunctie ongeveer lineair is. In dat geval is er een goede indicatie dat de benadering door de GPD geldig is. In het algemeen kan de gemiddelde overschrijdingsfunctie worden bestempeld als een goed hulpmiddel. Het hulpmiddel geeft echter geen unieke keus voor u.
Het POT-model heeft verder de volgende eigenschappen, zie [5]:
• Het aantal overschrijdingen volgt een Poisson-verdeling.
• Het maximum van de overschrijdingen is GEV-verdeeld.
Het POT-model heeft als belangrijke aanname dat de waarnemingen identiek en
onafhan-kelijk verdeeld zijn 1. Afhankelijkheden kunnen worden opgespoord door bijvoorbeeld een
scatterplot. Wanneer er afhankelijkheden zijn in de data kan er gekozen worden de waarne-mingen te groeperen, door bijvoorbeeld het maximum of het gemiddelde van een groep te nemen. De aanname dat de waarnemingen identiek verdeeld zijn kan worden gecontroleerd door middel van recordreeksanalyse ([6]):
Stel X = (X1, . . . , Xn)n∈N onafhankelijk en identiek verdeeld. Een waarneming Xi is een
record als geldt:
Xi > Mi−1.
Definieer nu het telprocess F door:
F1 = 1 Fn = 1 + n X j=2 I{Xj>Mj−1}. Er geldt nu: E[Fn] = n X j=1 1 j var(Fn) = n X j=1 1 j − 1 j 2 .
Het aantal records in de waarnemingen kan vergeleken worden met het aantal records wat verwacht wordt op basis van de recordreeksanalyse. Een grote afwijking is een indicatie dat de waarnemingen niet identiek verdeeld zijn.
Samenvattend: Onder de voorwaarden dat waarnemingen
1. in het maximum attractiedomein liggen van de GEV-verdeling,
2. onafhankelijk en identiek verdeeld zijn,
geldt voor alle waarnemingen die boven een bepaalde grens u komen:
1. dat de overschrijdingen boven u verdeeld zijn volgens de GPD
2. het maximum van deze overschrijdingen een GEV-verdeling heeft
3. dat deze overschrijdingen arriveren volgens een Poisson process.
Hoofdstuk 3: Schatters voor de GEV-verdeling en het
POT-model
In dit hoofdstuk zal ik kwantielschatters introduceren voor data die (bij benadering) af-komstig zijn van de GEV verdeling en voor data waarop het POT-model wordt toegepast. De kwantielschatters zijn een functie van de parameterschatters van de GPD-verdeling of de GEV-verdeling. In het algemeen kan worden gesteld dat het lastig is goede hoge kwantielschattingen te verkrijgen voor beide type verdelingen. De moeilijkheid wordt ver-oorzaakt door het beperkt aantal extreme waarnemingen. De onzekerheid in de schattingen is hierdoor groot.
Ik verdeel de kwantiel-schatters voor data (bij benadering) afkomstig van een GEV ver-deling in twee groepen. Er zijn schatters die veronderstellen dat de data afkomstig is van een GEV verdeling en er zijn schatters die veronderstellen dat de data afkomstig zijn van een verdeling uit het maximum attractie domein van de GEV verdeling. Dit laatste is een minder sterke veronderstelling dan de eerste veronderstelling. De tweede groep schatters kan dus vaker toegepast worden dan de eerste groep.
Tot de eerste groep schatters behoren de maximum-likelihood schatter (MLE) en een schat-ter gebaseerd op de probability weighted moments. Tot de tweede groep schatschat-ters behoren de Hill-schatter, Pickands-schatter en de Dekkers-Einmahl-De Haan schatter. De schatters in de tweede groep blijken minder makkelijk toe te passen dan die uit de eerste groep.
De eerste groep parameterschatters voor de GEV verdeling
Zoals ik in de inleiding opmerkte kunnen parameterschatters voor de GEV verdeling worden ingedeeld worden in twee groepen. Ik behandel eerst de schatters uit de eerste groep, ik maak daarom de volgende aanname:
Laat x = (X1, . . . , Xn) onafhankelijk en identiek verdeeld en zodanig dat geldt:
Xi ∼ G, i = 1, . . . , n,
waarbij G is de GEV verdeling.
De methode gebaseerd op de Maximum Likelihood
aan: l(ξ, ϕ, µ) = −n ln(ϕ) − (1 − ξ) n X i=1 yi+ n X i=1 expyi,
met yi = −ξ−1log(1 − ξ(xi− x0)/ϕ, zie ondermeer [12]. Voor het geval a priori de eis ξ = 0
wordt opgelegd reduceert de log-likelihood tot:
l(0, ϕ, µ) = −n ln(ϕ) n X i=1 exp −Xi− µ ϕ − n X i=1 −Xi− µ ϕ
Het opschrijven van de afgeleides naar de verschillende parameters en het gelijkstellen aan nul van deze afgeleides leidt niet tot een gesloten uitdrukking voor de parameters, daarom zijn er numerieke procedures ontwikkeld. Deze zijn inmiddels ge¨ımplementeerd in verschil-lende statistische programma’s, het numeriek bepalen van deze schatters levert geen grote problemen meer op. Aangezien in de meeste toepassingen aan voldaan is aan de eis dat ξ > −1/2, is maximum likelihood estimation een goede schatter voor de parameters van de GEV.
De methode van de probability-weighted moments
Een alternatieve schattingsmethode is de de methode van de probability-weighted mo-ments. Deze methode is een generalisatie van de klassieke momentenschatters. De methode is gebaseerd op het gelijkstellen van de theoretische momenten bepaald met behulp van de eigenschappen van de GEV aan de empirische momenten op basis van de data. Hosking, Wallis en Wood [13]laten in hun artikel het volgende zien:
Definieer
wr(θ) = E[XGrθ(X)], r ∈ N,
met Gθ de GEV verdeling en X een stochast met verdeling Gθ waarbij de parameter
θ = (ξ, µ, ϕ). De empirische tegenhanger van wr(θ) wordt gegeven door de vergelijking:
ˆ wr(θ) =
Z ∞
−∞
xGrθ(x)dFn(x), r ∈ N,
hierin is Fn de empirische verdelingsfunctie van de data X1, . . . , Xn. De basis voor de
Via ondermeer een kwantieltransformatie kan het volgende resultaat bewezen worden: wr(θ) = 1 n n X i=1 XjUjr,
waarin Xj is de j-de order statistic is en Uj de j-de order statistic uit een trekking uit de
uniforme verdeling ter grote van n. Er geldt: wθ = 1 1 + r µ −ϕ ξ(1 − Γ(1 − ξ)(1 + r) ξ , ξ 6= 0, ξ < 1.
Door het bepalen van de eerste 3 momenten kan de volgende vergelijking worden afgeleidt: 2w1(θ) − w0(θ)
3w2(θ) − w0(θ)
= 2
ξ− 1
3ξ− 1
Via het invullen van de empirische tegenhangers van de w0is aan de linkerzijde van de
vergelijking volgt een schatter voor ξ, deze noem ik ˆξ. Zodra deze bekend is kunnen de
andere schatters eenvoudig gevonden worden via de vergelijkingen: ˆ ϕ = (2 ˆw1− ˆw0) ˆξ Γ(1 − ˆξ)(2ξˆ− 1) ˆ µ = ˆw0+ ˆ ϕ ˆ ξ(1 − Γ(1 − ˆξ)).
Deze schattingsmethode is een ad-hoc methode, verder is er in de literatuur weinig be-kend over de theoretische eigenschappen van deze schatters. Computersimulaties geven aan dat voor data met minder dan 100 waarnemingen de probability weigthed moments vergelijkbaar zijn in effici¨enty met de maximum likelihoodschatters. Wanneer het aantal waarnemingen toeneemt is MLE effici¨enter. De bias van de probability-weigthed moments blijkt klein in computersimulaties en daalt snel naarmate het aantal waarnemingen toe-neemt. Een ander groot voordeel van deze schatters boven andere schatters is dat ze snel te bepalen zijn en bovendien altijd nette uitkomsten geven. Dat wil zeggen, er worden geen schattingen gevonden die onzinnig zijn, bijvoorbeeld ϕ < 0. Nadat er schattingen zijn verkregen voor de parameters van de GEV verdeling is het relatief eenvoudig een
kwantiel-schatter te vinden. Laat xp het p-de kwantiel van de GEV-verdeling, oftewel:
xp := G−1(p),
met G−1(p) = µ − ϕξ(1 − (− ln(p)−ξ), de inverse van de GEV verdeling. Een natuurlijk
schatter ˆxp voor xp is nu:
ˆ xp = ˆµ − ˆ ϕ ˆ ξ h 1 − (− ln(p)− ˆξ)i.
De tweede groep parameterschatters voor de GEV verdeling
De tweede groep kwantiel-schatters voor de GEV verdeling is ontwikkeld voor data x waarvoor het volgende geldt:
(X, X1, . . . Xn) = x ∈ M DA(G),
waarbij G is de GEV verdeling. Het kan bewezen worden (zie Goldie en Resnick [15]) dat de volgende benadering werkt voor een grote u:
n ¯G(u) ≈ 1 + ξu − dn cn −1ξ , (2)
waarbij cn en dn de normalizerende constantes zijn uit hoofdstuk 1 zoals in vergelijking 1.
Voor een kwantiel-schatter ˆxp gebaseerd op bovenstaande benadering zijn nu schatters voor
cn, dn en ξ noodzakelijk. Voor ξ bestaan Pickands-schatter, Hills-schatter en de schatter
van Dekkers, Einmahl en De Haan. Al deze schatters zijn gebaseerd op ’order statistics’ uit de data. Om deze schatters toe te passen moet een keuze gemaakt worden over het aantal order statistics dat gebruikt zal worden. Als hulpmiddel bij deze keuze zijn verschillende instrumenten beschikbaar (bijvoorbeeld de Hill-plot), deze instrumenten geven echter geen direct sluitend antwoord voor de te maken keuze. De keuze die gemaakt moet worden wordt een moeilijke door het volgende dilemma.
De zojuist ge¨ıntroduceerde benadering is alleen geldig in de staart van de verdeling. Wan-neer teveel order statistics gebruikt worden, is er een kans dat de benadering niet meer geldig is en een bias in de schattingen zal sluipen. Als echter te weinig data gebruikt worden, neemt de variantie in de schatters toe. Ik geef in deze paragraaf de verschillende schatters voor ξ en enkele van hun eigenschappen.
Pickands Schatter
Pickands schatter voor ξ kan worden bepaald met de volgende vergelijking: ˆ ξ(n−k+1)P = 1 ln 2ln X(n−k+1)− X(n−2k+1) X(n−2k+1)− X(n−4k+1) ,
waarbij Xj is de j-de order statistic. In deze schatter moet een keuze voor k gemaakt
De eigenschappen van Pickands schatter (Dekkers en De Haan[17]) zijn als volgt:
1. Zwak consistent, als k → ∞, k/n → 0 terwijl n → ∞ dan convergeert ˆξP in
waar-schijnlijkheid naar ξ.
2. Sterk consistent, als k/n → 0, k/ ln(ln(n)) → ∞ terwijl n → ∞ dan convergeert ξ
met kans 1 naar ξ.
3. Onder ingewikkelder technische voorwaarden voor k en de verdelingsfunctie F van
de data geldt: p(k)(ˆξ − ξ)→ N (0, v(ξ)), als n → ∞, met v(ξ) =d (2(2ξ2(2ξ−1) ln(2))2ξ+1+1)2. Deze eigenschappen betekenen in de praktijk dat k groot gekozen moet worden, maar ook zo dat n − k groot blijft. Deze ’sloppy’ formulering geeft aan dat er nogal wat haken en ogen aan de keus voor k zitten. De eigenschappen vormen dan ook zeker geen garantie voor een betrouwbare schatter, de keuze van k is een moeilijke en de vraag is of er voldoende waarnemingen bestaan om een betrouwbare schatter te verkrijgen. ’Pickands plot’ kan een
handig hulpmiddel zijn. Dit is een plot waarbij Pickands schatter ˆξ wordt uitgezet tegen
k. Voor elke waarde van k wordt dan Pickands schatter ˆξ uitgerekend. Het is te hopen dat
naarmate k toeneemt, de waarden van ˆξ stabiliseren. De waarde van k moet niet gekozen
worden uit een interval waarop de Pickands schatter ˆξ sterk varieert. Hangt de waarde van
Pickands schatter sterk af van de waarde k, dan is dit een indicatie dat k te groot is.
Hills schatter
Op verschillende manieren (maximum likelihood, mean excess functie) is een afleiding te vinden voor Hills schatter, deze is ontworpen voor ξ > 0 en ongeschikt voor andere waarden van ξ, zie Hill [16]. Definieer α door α := 1
ξ, Hills schatter heeft dan de volgende vorm:
ˆ αHk = 1 k k X j=1 ln(Xn−j+1) − ln(Xn−k) !−1
Onder bepaalde voorwaarden gelden voor Hills schatter dezelfde eigenschappen als die voor Pickands’, zowel zwak als sterk consistent en assymptotisch normaal verdeeld. Net als bij Pickands schatter is de keuze voor k een lastige. Hier moet de eerder genoemde afweging gemaakt worden tussen het aantal ’order statistics’ en de geldigheid van benadering 2. Een Hill-plotkan een handig hulpmiddel zijn; een evaluatie van ˆαH
k voor verschillende waarden
Dekkers, Einmahl en De Haan
Een andere schatter voor ξ in deze context is die van Dekkers, Einmahl en De Haan. Dit is een uitbreiding van Hills schatter in die zin dat deze schatter geschikt is voor alle ξ terwijl Hills schatter slechts geschikt is voor ξ > 0. De schatter is gedefinieerd door:
ξkDEH = 1 + Hn1 +1 2 (H1 n)2 Hn(2) − 1 −1 , waarbij H1 n = 1 k Pk j=1ln(Xn−j+1) − ln(Xn−k) en H2 n = 1 k Pk j=1ln(Xn−j+1) − ln(Xn−k) 2 .
Deze schatter heeft vergelijkbare eigenschappen als de vorige twee schatters.
Nadat er drie schatters zijn gepresenteerd dringt zich de vraag op: ’Welke is de beste?’. Net als met de keuze voor de hoeveelheid order statistics die gebruikt moeten worden bij het toepassen van de schatters, is ook op deze vraag geen pasklaar antwoord te geven. In het algemeen kan worden gezegd dat Hills schatter ongeschikt is voor ξ ≤ 0, dat Pickands schatter afhankelijker is van de keuze voor het aantal order statistics dan de andere twee. Ik naast de twee schatters uit de vorige paragraaf (ML-schatter, pwm-schatter) voornamelijk gebruik maken van Pickands schatter omdat deze relatief gemakkelijk te implementeren is.
De normaliserende constanten
We hebben nu 3 schatters voor de parameter ξ. Voor een kwantiel-schatter zijn echter ook nog schatters nodig voor de normaliserende constanten. De normaliserende constan-ten hangen in het geval van verdelingen die in het maximum attractie domein van de GEV verdeling liggen af van de kwantielen in de staarten van deze verdelingen. Dit ter-wijl deze kwantielen juist geschat moeten worden. Dit maakt het bijzonder moeilijk een kwantiel-schatter te vinden. Voor het geval ξ > 0 is er echter een methode. Via een slimme transformatie (zie Dekkers en Haan [17]) kan echter de volgende kwantiel-schatter worden gevonden: ˆ xp = n k(1 − p) − ˆξH n,k ∗ Xn−k
Balans
Ik heb nu een kwantiel-schatter ge¨ıntroduceerd voor het geval dat de data afkomstig zijn van de GEV verdeling. Deze veronderstelling over de data is in de praktijk moeilijk te verifi¨eren. Het is niet bekend of de data netjes verdeeld zijn volgens de GEV. De betrouwbaarheid van deze schatter is dan ook lastig in te schatten. Ondermeer vanwege dit bezwaar zijn er andere schatters ontwikkeld. De veronderstelling dat de data komen van een verdeling die in het maximum attractiedomein van de GEV ligt zal vaker geldig zijn. Voor dit geval heb ik drie schatters ge¨ıntroduceerd voor ξ en een kwantiel-schatter die echter beperkt geldig is.
Schatters voor het POT-model
In hoofdstuk 1 heb ik het POT-model ge¨ıntroduceerd, in deze paragraaf zal ik beschrijven hoe kwantielen in dit model kunnen worden geschat. Dit model is geschikt voor data gege-nereerd uit een verdeling die ligt in het maximum attractiedomein van de GEV verdeling. Met andere woorden:
(X1, . . . , Xn) := x ∈ M DA(G), waarbij G is de GEV verdeling.
Ik definieer vervolgens: ¯F (x) = 1 − F (x). De basis voor de maximum likehoodschatter en
de probability weighted moments is de volgende vergelijking: ¯
F (u + y) = ¯F (u) ¯Fu(y),
waarbij ¯Fu(y) = P (X − u > y|X > u).
Definieer nu Nu :=Pni=1IXi>u. Een natuurlijk schatter voor ¯F (u) wordt gegeven door:
ˆ ¯
F (u) = Nu/n.
Laat Pξ,β(x) een GPD met parameters β en ξ. Vanwege het feit dat overschrijdingen over
een hoge grens GPD verdeeld zijn, wordt een schatter voor ¯F (u + y) nu gegeven door:
ˆ¯
Fu(y) = Pξ, ˆˆβ(x).
De dichtheidsfunctie behorende bij de GPD is: p(x) = 1 β 1 + ξx β −1ξ−1 , x ∈ D(ξ, β), met D(ξ, β) =( [0, ∞)h als ξ ≥ 0 0,−βξ i als ξ < 0
Stel x = (X1, . . . , Xn) waarbij geldt: Xi ∼ GP D, met de Xi’s onafhankelijk van elkaar. De
log-likelihood l bij de data x als functie van de parameters ξ en β is: l((ξ, β); x) = −n ln(β) − 1 ξ + 1 n X i+1 ln 1 + ξ βXi I{Xi∈D(ξ,β)}.
Deze likelihoodfunctie kan geoptimaliseerd worden via numerieke methoden. Omdat het domein van de likelihood afhankelijk is van de parameters is voorzichtigheid geboden. Smith [18] laat zien dat de mle goed werkt voor ξ > −1/2. In dat geval geldt:
n1/2 ξˆn− ξ, ˆ βn β − 1 ! d → N (0, M−1), voor n → ∞ met M−1 = (1 + ξ) 1 + ξ −1 −1 2 .
De methode van probability weighted moments (Hosking [14]) maakt gebruik van de vari-abelen:
wr = E[Z( ¯Gξ,β(Z))r], voor r = 0, 1.
Dit leidt direct tot:
β = 2w0w1
w0− 2w1
en ξ = 2 − w0
w0− 2w1
.
De theoretische variabelen w0 en w1 vervangen door hun empirische tegenhangers geeft
ad-hoc schatters voor β en ξ. De schatters staan in de literatuur bekend als een redelijk alternatief voor de maximum likelihood methode voor het typische geval dat: ξ ≥ 0. De betrouwbaarheid van de schatters is afhankelijk van de verdeling F maar ook van de keuze voor u. Zoals in hoofdstuk 1 reeds opgemerkt is hierbij een handig middel de mean-excess functie. De kwantiel-schatter kan nu gemakkelijk worden verkregen door het invullen van de schatters in de vergelijking: ˆ xp = u + ˆ β ˆ ξ n Nu (1 − p) − ˆξ − 1 ! .
Hoofdstuk 4: Copula’s
In dit hoofdstuk zal ik algemene theorie over copula’s beschrijven. Allereerst introduceer ik formeel een copula en demonstreer ik hoe deze gebruikt kan worden bij het beschijven twee-dimensionale verdelingen. Daarna beschrijf ik Kendalls correlatie en introduceer ik de Archimedische familie. Tot slot geef ik nog enkele voorbeelden van dit type copulafuncties.
Introductie copula’s
Een copula kan gebruikt worden om een gezamenlijke verdelingsfunctie van twee of meer variabelen weer te geven. Een copula verdeelt de gezamenlijke verdelingsfunctie in twee bijdrages, een bijdrage van twee marginale verdelingen en een bijdrage van de afhanke-lijkhheidsrelatie tussen de twee variabelen. Venter [7] geeft de volgende definitie voor een copulafunctie:
Neem twee stochasten X en Y met corresponderende verdelingsfuncties FX(x) en FY(y). De
gezamenlijke verdelingsfunctie FX,Y(x, y) kan op de volgende manier worden weergegeven
met behulp van een copula C:
FX,Y(x, y) := C(FX(x), FY(y)).
Definieer C(i, j) door C(i, j) = FX,Y(FX−1(i), F −1
Y (j)). Deze definitie betekent dat de
vol-gende vergelijking opgesteld kan worden:
C(FX(x), FY(y)) = FX,Y(FX−1(FX(x)), FY−1(FY(y))) = FX,Y(x, y),
oftewel C(FX(x), FY(y)) = P (X ≤ x, Y ≤ y). Met andere woorden, C(FX(x), FY(y)) is
gelijk aan de gezamenlijke verdelingsfunctie van de stochasten X en Y. Merk op dat het domein van de copulafunctie gelijk is aan: I2 = [0, 1]X[0, 1], en het bereik gelijk is aan:d
I1 d= [0, 1]. Een voorbeeld van een simpele copulafunctie is:
C(i, j) = ij.
Deze copulafunctie geeft de gezamenlijke verdeling van twee onafhankelijke stochasten.
Twee andere simpele copulafuncties zijn de zogenaamde Frech´etlimieten, zie ondermeer
Wang [8]:
M (i, j) := min(i, j)
W (i, j) := max(i + j − 1, 0) Er geldt:
W (i, j) ≤ C(i, j) ≤ M (i, j).
Kendalls correlatie
Copulafuncties zijn in staat afhankelijkheden tussen stochastische variabelen te beschrij-ven. In tegenstelling tot de bekende standaardmaat ’(lineaire) correlatie’, die een maatstaf is voor de algemene afhankelijkheid tussen twee stochasten, bieden copulafuncties de mo-gelijkheid tot een dynamische beschrijving van de afhankelijkheidsrelatie. Dit kan nuttig zijn in tal van voorbeelden. De relatie tussen twee variabelen in de staarten van de verde-lingen kan verschillen met de relatie tussen dezelfde variabelen rond hun verwachtingen. De afhankelijkheid tussen bijvoorbeeld hoeveelheid neerslag en de opbrengst per vierkante meter landbouwgrond is waarschijnlijk in de staarten van de verdeling sterker dan rond het centrum van de verdeling. Zodra er extreme hoeveelheden neerslag vallen zal de opbrengst waarschijnlijk extreem laag zijn. Een hele droge of hele natte periode is in het algemeen desastreus voor de opbrengst. Er is in de staart dus een sterke relatie tussen beide groot-heden.
Vallen er neerslaghoeveelheden die niet teveel afwijken van het gemiddelde dan is het waar-schijnlijk dat de opbrengst per vierkante meter landbouwgrond minder afhankelijk is van de neerslag, maar meer van andere factoren. In het centrum van de verdeling is de corre-latie dan veel minder sterk. Een copulafunctie geeft de mogelijkheid om deze dynamische afhankelijkheid te modelleren.
Een interessant punt om op te merken is het volgende: twee verschillende paren stochas-ten kunnen dezelfde copulastructuur hebben maar niet dezelfde (lineaire) correlatie. Een alternatieve correlatiemaat die wel constant is over copulafuncties is de Kendall correlatie. Deze maat zal in hoofdstuk 5 een belangrijke rol spelen in de schattingsprocedure voor de parameters in de te gebruiken Copulafunctie. Er zijn verschillende definities voor de Kendall correlatie (ρ), zie Nelson, [10], ik geef een simpele uit Ventor:
ρ = 4E[C(i, j)] − 1 = 4
Z Z
I2
C(i, j)dC(i, j)) − 1. Voor de simpele copulafuncties van perfect gecorreleerde variabelen geldt:
C(i, j) = i = j, dus E[C(i, j)] = Z 1 0 Z 1 0 i di dj = 0.5,
waaruit volgt dat ρ = 1. Voor onafhankelijke variabelen daarentegen geldt: C(i, j) = ij, dus E[C(i, j)] = Z 1 0 Z 1 0 i di dj = 1/4,
waaruit volgt dat ρ = 0. Het kan bewezen worden met behulp van de Frech´et-limieten voor
De schattingsprocedure die ik ga gebruiken om een copulafunctie te schatten maakt ge-bruik van de zogeheten Archimedische copulafamilies. Een copulafamilie bestaat uit alle copulafuncties die slechts in parameterwaarde(n) verschillen. Voor het introduceren van Archimedische copulafamilies maak ik gebruik van de definitie van Cherubini [9]:
Laat φ een functie van [0, 1] naar [0, ∞] die convex, continu en strict dalend is en waarvoor bovendien geldt: φ(1) = 0. De pseudo-inverse van φ is nu gedefinieerd door:
φ−1(t)=d φ
−1(t) 0 ≤ t ≤ φ(0)
0 φ(0) ≤ t ≤ ∞
Een copulafamilie C(i, j) is van het Archimedische type als geldt:
C(i, j) = φ−1(φ(i) + φ(j)) (3)
Bekende copulafamilies die tot deze familie behoren zijn Gumbels copulafamilie en de Claytons copulafamilie (ook wel Heavy Right Tail genoemd), zie o.a. Genest en Rivest. Hieronder geef ik de definities van deze copulafamilies 2:
Gumbel: CG(i, j) = exp
n
−[(− ln(i))a+ (− ln(j))a]1ao, a > 0.
Clayton: CCl(i, j) = [i−a+ ja− 2]−
1
a, a > 0.
Tot slot formuleer ik nog een belangrijk resultaat uit Genest en Rivest:
Veronderstel CA(i, j) een Archimedische copulafunctie en veronderstel X en Y zijn uniform
verdeelde variabelen. Definieer nu λ(v) = φφ(v)0(v). Er geldt nu:
De grootheid V = CA(X, Y ) volgt de verdeling K(v) := v − λ(v)
De verdeling K(v) staat bekend als ’Kendall’s process’ en speelt een belangrijke rol in het volgende hoofdstuk.
2. Voor de Gumbel verdeling geldt bijvoorbeeld: φ(i) = (− ln(i))a en dus φ−1(t) = exp(−t1
a), voor een
Hoofdstuk 5: Modelleren met Copula’s
In dit hoofdstuk zal ik een beperkt overzicht geven van de huidige inzichten in het model-leren met Copula’s. Er zijn veel papers geschreven over dit onderwerp en het laatste woord hierover is nog lang niet gesproken. Een belangrijk probleem bij het gebruik van copula’s is het selecteren van een geschikte copulafamilie. Er is nog geen unieke procedure bekend welke altijd goede resultaten oplevert.
Durrleman, Nikeghbali en Roncalli [19] onderscheiden drie verschillende situaties waarin copula’s bij het modelleren gebruikt worden. De eerste situatie gaat uit van de veronder-stelling dat (bijvoorbeeld op theoretische gronden) de copulafamilie bekend is, maar dat de parameter(s) nog bepaald moet worden. De tweede situatie is niet gebaseerd op een para-metrische copulafunctie, maar gebaseerd op een onbekende onderliggende afhankelijkheids-structuur. De derde situatie veronderstelt dat de copulafunctie geselecteerd moet worden uit een aantal copulafamilies. Met name de derde situatie zal ik uitgebreid bespreken, want deze zal ik toepassen in mijn onderzoek.
Een bekende copulafunctie
De eerste situatie veronderstelt het volgende:
Er zijn waarnemingen x1, . . . , xnen y1, . . . , yn. Deze zijn afkomstig van een marginale
verde-lingen FX(x, θx) en FY(y, θy). Deze verdelingen zijn afkomstig van een bekend
veronderstel-de familie waarvan alleen veronderstel-de bijbehorenveronderstel-de parameters θx en θy onbekend zijn. Daarnaast is
de afhankelijkheid tussen de waarnemingen te beschrijven met een bekend veronderstelde copulafamilie C(FX(x, θx), FY(x, θy), θc) waarbij θc is de vector met de onbekende
para-meter(s) van de copulafunctie. Er zijn twee methodes om de onbekende parametervector θ = (θx, θy, θc) te schatten. De eerste is gebaseerd op de log-likelihood waarbij alle
onbe-kende parameters simultaan geschat worden. Het intensieve en complexe rekenwerk van deze schattingsmethode kan een probleem opleveren. De tweede methode is ook gebaseerd op de log-likelihood, maar deze schat eerst de parameters van de marginale verdelingen en vervolgens de parameter(s) van de copulafunctie gebruik makend van de parameter-schattingen van de marginale verdelingen. Deze methode staat bekend als de IFM, een uitgebreide beschrijving is te vinden in Joe en XU [20].
Niet-parametisch schatten
Deze methode veronderstelt zoals eerder opgemerkt een onderliggende afhankelijkheids-structuur tussen twee marginale verdelingen. Deze afhankelijkheids-structuur is uniek. Om de afhankelijkheids-structuur te benaderen wordt een empirische copulafunctie bepaald. Deze manier van afhankelijkheden modelleren is consistent, over de kwaliteit van de eigenschappen is verder niet veel bekend. Ik zal deze methode niet toepassen op data.
Copula selectie
In deze paragraaf zal ik beschrijven hoe een keuze kan worden gemaakt tussen verschillende copulafamilies. In de literatuur is veel aandacht voor de reeds ge¨ıntroduceerde Achimedische copulafamilies. Waarschijnlijk is dit het geval omdat er relatief eenvoudige procedures bestaan om voor dit geval een copulafunctie te bepalen. Ik sluit me bij deze praktijk aan en beperk mij tot deze copulafamilies. Genest en Rivest demonstreren [22] een methode waarmee een keuze gemaakt kan worden tussen Archimedische copulafamilies. Ik behandel in deze paragraaf de methode voor een twee-dimensionale copulafunctie. De methode kan worden uitgebreid naar hogere dimensies. Aan het einde van dit hoofdstuk geef ik enkele suggesties om deze procedure te verbeteren.
Veronderstel we hebben waarnemingen (X1, Y1), (X2, Y2), . . . , (Xn, Yn) waarvoor geldt:
Xi ∼ FX, i = 1, . . . , n
Yi ∼ FY, i = 1, . . . , n.
Uit hoofdstuk 4 herhaal ik Kendalls process, de univariate verdelingsfunctie K(v):
K(v) = P (H(X, Y ) < v) = P r(C(F (X), F (Y )) < v) = v − λ(v), (4)
waarin λ(v) = φφ(v)0(v). De schatter die nu volgt maakt gebruik van het feit dat de empiri-sche tegenhanger, de grootheid Hn(Xi, Yi), het percentage waarnemingen (Xj, Yj) bevat
waarvoor geldt dat:
Met behulp van dit feit kan de empirische verdelingsfunctie van Kendalls proces uitgerekend
worden worden en aan de hand van de volgende vergelijking3 :
ˆ K(v) = 1 n n X i=1 I[θ1i<v] , (5) met θi = 1 n − 1 n X j=1
I[X(j)<X(i),Y (j)<Y (i)]2
en met I1
x<a een indicatorfunctie gedefinieerd door:
I[x<a]1 = 1 als x < a
0 elders ,
en een analoge definitie voor indicatorfunctie I2.
Met behulp van vergelijking 4 voor K(v) kan een uitdrukking voor λ(v) bepaald worden. Via de differentiaalvergelijking voor λ(v) kan ook φ(v) bepaald worden. Omdat φ de ka-rakteriserende eigenschap is van iedere Archimedische copulafamilie kan de juiste familie worden geselecteerd4.
Een iets slimmere methode is theoretische functie K(v) behorende bij verschillende copula-functies te bepalen op grond van de vergelijking K(v) = v − λ(v). Vervolgens kan deze
vergeleken worden met de empirische waarde, Kn(v). Dit bespaart enig rekenwerk omdat
de functie φ(v) niet geheel bepaald hoeft te worden. Een overzicht van de relevante theo-retische grootheden behorende bij de verschillende copulafamilies staat in tabel 1.
Op dit moment heb ik de strategie beschreven om tot een selectie van de copulafamilie
te komen. De strategie brengt ´e´en probleem met zich mee. Er moet namelijk een keus
gemaakt worden voor de parameter(s) van de copulafunctie. Zoals blijkt uit tabel 1 heeft iedere copulafamilie minimaal 1 parameter, maar er zijn ook copulafuncties die afhangen van meer parameters zoals copulafuncties uit de log-copulafamilie. Om de empirische
waar-den van Kendalls process Kn(v) te vergelijken met die van de theoretische waarden K(v)
zal allereerst een waarde voor deze parameters gekozen moeten worden. E´en manier om dit
te doen is door voor de verschillende copulafamilies een maximum likelihood-schatting te doen. Deze methode heb ik globaal omschreven in dit hoofdstuk. Een andere aanpak die
werkt in het geval slechts ´e´en parameter gekozen dient te worden is de momentenschatter.
Dit houdt in dat de parameter zo gekozen wordt dat Kendalls correlatie uit de steekproef gelijk is aan de theoretische waarde. De theoretische waarden staan in tabel 1. Deze mo-mentenschatter is zeer simpel voor de Clayton- en Gumbelcopulafamilie, voor de Frank-en Log-copulafamilie betekFrank-ent dit eFrank-en lastige numeriek procedure. GFrank-enest Frank-en Rivest [22] kiezen voor deze laatste momentenschatter, ik zelf suggereer later in dit hoofdstuk nog een
3. Merk op het subtiele verschil tussen θ(i) en Hn(Xi, Yi). Genest en Rivest merken op dat θ geschikter
is dan Hn(Xi, Yi) omdat het bereik van de eerste grootheid gelijk is aan het interval (0,1).
4. φφ(v)0(v)= λ1 ⇒ φ(v) = exp(Rv
v0
andere procedure.
Kendalls correlatie kan uit de steekproef worden gehaald door middel van de vergelijking: τn= 4 ¯V − 1,
met ¯V =Pn
i=1θi
Tabel 1: Overzicht eigenschappen copulafamilies .
Copulafamilie φ(v) λ(v) τ
Clayton (v−a− 1)/a v(1 − va)/a a/(a + 2)
Frank log 1−exp(−α) 1−exp(−αv) 1−exp(−αv) α exp(−αv) log( 1−exp(−α) 1−exp(−αv) 1 + 4(D1(a)−1) a
Gumbel (− log(v))α+1) −v log(v)/(α + 1) a/(a + 1)
Log-copula (1 −log(v)ag )a+1− 1 agv[(1−log(v)/ag)a+1−1]
(a+1)(1−log(v)/ag)a
a−2+4C(a) (a+1)g
D1is de Debye function van order 1,R a 0(t/a(e t−1))dt en C is de exponenti¨ele integraalR∞ 0 (e −at/(1+t)a)dt, ik ga niet in op de details.
De procedure om een keus te maken tussen verschillende Archimedische families is samen-gevat als volgt:
1. Bereken Kendalls tau op basis van de data.
2. Bepaal parameterschatters voor iedere familie op grond van Kendalls tau.
3. Bereken het empirische Kendalls proces Kn(v).
4. Vergelijk Kn(v) met de theoretische functies K(v) en kies de copulafamilie waarvan
Kendalls proces K(v) het beste Kn(v) benaderd.
Genest en Rivest kiezen ervoor de laatste stap uit te voeren door middel van visuele in-spectie, Frees en Valdez [21] kiezen voor een vergelijking op basis van een QQ-plot of voor een iets gewijzigde procedure. Deze is als volgt:
1. Bereken de empirische verdelingsfunctie Kn(v).
2. Minimaliseer over alle families en mogelijke parameterwaarden de afstand tussen
Het minimaliseren van de afstand kan op basis van de discrete L2 norm. Definieer Kθii(v)
de theoretische tegenhanger van de empirische functie Kn(v) behorende bij copulafamilie
i en parametervector θi. Definieer Θj als de parameterruimte waarin θj zich kan begeven.
Formeel ziet deze minimalisatie er zo uit: min j=1,...kθminj∈Θj t X i=1 Kn(i/t) − Kθjj(i/t) 2
In deze formule is t het aantal discrete punten waarop de theoretische en empirische functie vergeleken worden.
Goodness of Fit
Een interessante vraag is vervolgens hoe goed deze modellen zijn. Hierbij zijn twee midde-len die kunnen helpen. De eerste is een Pearson chi-squared goodness-of-fit (GOF) toets, waarbij waargenomen frequenties vergeleken worden met de frequenties op basis van het model. De toets werkt als volgt:
Eerst maak ik groepen met als grenzen de ’order statistics’ (X([i∗n/k]), Y([j∗n/k])), i =
1, . . . , k; j = 1, . . . , k met k de wortel van het aantal groepen, n het aantal waarnemingen en waarin [m] is de afgeronde waarde van m. Definieer de waargenomen waarde O(i, j) voor i, j = 1, . . . , k door:
O(i, j) := card(Xi, Yi) : X([(i−1)∗n/k])< Xi < X([i∗n/k]), Y([(i−1)∗n/k]) < Yi < Y([i∗n/k]) ,
Met card de cardinaliteit van de verzameling. Definieer de verwachte waarde E(i, j) door: E(i, j) := Cα(i/k, j/k) ∗ n; i, j = 1, . . . , k,
waarbij Cα(i, j) is de geschatte copulafunctie ge¨evalueerd in (i,j). De uitdrukking voor
deze copulafunctie wordt verkregen door de karakteristieke functie φ(v) uit tabel 1 te substitueren in vergelijking 3. De toetsgrootheid
X := k X i=1 n X j=1 (O(i, j) − E(i, j))2 E(i, j) ,
De variantie van Kendalls proces K(v), nodig voor de constructie van betrouwbaarheids-intervallen, kan verkregen worden via een gecompliceerde afleiding en levert een
gecompli-ceerde uitdrukking op. De variantie van Kendalls proces noem ik voortaan σ2
α(v). Alleen
voor de familie van Clayton bestaat een gesloten uitdrukking voor de variantie van K(v). Deze gesloten uitdrukking wordt door Genest en Rivest ook gebruikt als geschikte bena-dering voor de variantie van Kendalls proces van andere Archimedische families. Voor de afleiding van de gesloten uitdrukking verwijs ik naar de appendix van Genest en Rivest [22], ik geef hier alleen de formule:
σα2 = K(v) ¯K(v) + k(v)k(v)R(v) − 2v ¯K(v) /n , met
¯
K(v) = 1 − K(v), k(v) = K0(v) en
R(v) = EC {min(X1, X2), min(Y1, Y 2)} − v2|C(X1, Y1) = C(X2, Y 2) = v .
Voor de copulafamilie van het type Clayton zijn de afzonderlijke componenten in gesloten uitdrukkingen op te schrijven: K(v) = v(1 + (1 − vα)/α) k(v) = (α + 1)(1 − vα)/α R(v) = 2αv (1 − α)(1 − 2α)(1 − vα)2 [α(2 − vα)(1 − 2α) − α] − v2.
Het 95%-betrouwbaarheidsinterval wordt nu gegeven door:
[Kn(v) − 4.72 ∗ σ2α, Kn(v) + 4.72 ∗ σα2] (6)
Hoofdstuk 6: Data en Resultaten
Inleiding
In dit hoofdstuk presenteer ik de resultaten van mijn analyse. De dataset die ik onder-zocht heb bestaat uit dagelijkse waarnemingen met betrekking tot zonneschijn, neerslag, temperatuur,windsnelheid en luchtdruk gedaan op het weerstation ’De Bilt’. De dataset is geheel belangeloos beschikbaar gesteld door het ’KONINKLIJK NEDERLANDS METE-OROLOGISCH INSTITUUT’(KNMI) en is uiteraard verkrijgbaar bij de auteur van dit stuk, maar ook te vinden op de website van het KNMI. De gebruikte waarnemingen lopen vanaf 1 januari 1951 tot en met 31 december 2006. Waarnemingen van eerdere tijdstippen zijn wel beschikbaar voor enkele type waarnemingen, maar niet voor alle typen.
De resultaten van de analyse bestaan uit twee groepen. De eerste groep resultaten die ik presenteer is afkomstig van de analyse op basis van de extreme waarden theorie voor iedere afzonderlijke reeks waarnemingen. De tweede groep gepresenteerde resultaten is het resul-taat van het toepassen van de theorie over copula’s op de modellen voor de afzonderlijke reeksen. Bij deze analyse is het van belang dat iedere reeks even lang is, met het oog hierop heb ik ook iedere reeks evenlang gekozen.
Resultaten toepassing extreme waarden theorie
Allereerst presenteer ik de de resulaten van de analyse naar aanleiding van de extreme waarden theorie uit hoofdstuk 2 en 3. De theorie heb ik toegepast op 5 reeksen, bestaande uit de volgende type dagelijkse waarnemingen:
• De som van de neerslag in 0,1 millimeter.
• Het aantal uren zon in 0,1 uur.
• De windsnelheid in 0,1 m/s.
• De gemiddelde luchtdruk in 0,1 millibar.
• De gemiddelde temperatuur in 0,1 graden celsius.
Iedere reeks waarnemingen heb ik achtereenvolgens onderzocht op de volgende punten:
• Records.
• Onafhankelijkheid tussen waarnemingen.
• Verdeling van de extremen.
Voor het eerste punt heb ik gebruikt gemaakt van de recordanalyse uit het vorige hoofdstuk, voor het tweede punt heb ik gebruik gemaakt van visuele inspectie van een scatterplot. Wanneer de data onafhankelijk en stabiel verdeeld zijn kunnen de extremen onderzocht worden. Ik gebruik hiervoor de twee verschillende methodes uit het hoofdstuk over extre-me waarden theorie.
Dit doe ik enerzijds vanwege het feit dat deze gemakkelijker te gebruiken en te implemen-teren zijn, anderzijds omdat de kwaliteit van andere schatters zoals Hill en de schatter van Dekkers, Einmahl en De Haan afhangt van een arbitraire keus voor het aantal ’order statistics’. De maximum likelihoodsschatter en de methode van de probability weighted moments blijken goede resultaten te geven, dit meet ik af aan de mate waarin de theo-retische verdeling een goede beschrijving geeft van de empirische verdeling. De keus voor de makkelijk toepasbare schatters ligt hiermee voor de hand. De tweede methode is ge-baseerd op de theorie dat de data die afkomstig zijn van een verdeling uit het maximum attractiedomein van de GEV-verdeling boven een bepaalde grens GPD-verdeeld zijn. Dit zijn de zogenaamde POT-modellen waarvan ik de eigenschappen in het vorige hoofdstuk heb genoemd. Voor iedere reeks zal ik een mean-excess plot bepalen waarmee een grens u gekozen kan worden. De parameters van de GPD-verdeling schat ik vervolgens met de eerder ge¨ıntroduceerde schatters voor het POT-model, deze geven over het algemeen goede modellen.
In dit hoofdstuk bespreek ik de belangrijkste resultaten voor iedere reeks. Alle figuren waarnaar ik in dit hoofdstuk verwijs zijn te vinden in de bijlage.
De som van de hoeveelheid neerslag per dag
Een belangrijke eigenschap van de reeks waarnemingen van de hoeveelheid neerslag op een dag is dat ze allemaal groter of gelijk zijn aan nul. Een andere eigenschap van deze reeks is dat een groot percentage van de waarnemingen in deze catagorie ligt. Dit ligt voor de hand, er kunnen geen negatieve neerslaghoeveelheden vallen en in Nederland is het aantal droge dagen groter dan het aantal natte dagen.
Daarnaast valt op dat hoe groter de hoeveelheid neerslag, hoe minder vaak deze hoeveel-heid gesignaleerd wordt. De vorm van een histogram van deze reeks lijkt er op te wijzen dat de som van de neerslag per dag voor iedere dag een verdeling volgt met een eindig rechtereindpunt of eventueel een verdeling met ’normale’ staarten. De veronderstelling dat de neerslaghoeveelheid per dag identiek verdeeld is over de gehele periode gecombineerd met het aantal waarnemingen van ruim 3700 betekent dat de verwachting van het aantal records gelijk is aan 11,1 en de variantie van het aantal records gelijk aan 9,5. Dit resultaat levert geen bewijs op tegen de veronderstelling dat er de neerslaghoeveelheden identiek verdeeld zijn oftewel dat iedere dag dezelfde kansverdeling heeft met betrekking tot de hoeveelheid neerslag die zal vallen.
in onderstaande tabel. In de eerste kolom staat de betreffende reeks, in de resp. tweede en derde kolom het verwachte aantal en de variantie van het aantal records op grond van het aantal waarnemingen. In de vierde kolom staat het waargenomen aantal records. Van klimaatverandering op het gebied van neerslag blijkt uit deze analyse niets.
Reeks Verwacht Variantie Waargenomen
Som dagelijkse neerslag 10,5 8,8 9
Jaarlijks maximum 4,6 3,0 4
Een simpele scatterplot geeft geen aanleiding te veronderstellen dat er sprake is van afhan-kelijkheden zodat de extremen van deze reeks onderzocht kunnen worden met behulp van de extreme waarden theorie.
Allereerst modelleer ik de ruwe data aan de hand van eerder genoemde methoden. De eer-ste methode, het schatten van de parameters van de GEV-verdeling, geeft voor de reeks ruwe neerslaghoeveelheden een aantal problemen. In de eerste plaats zijn de empirische verdeling en de verdeling op basis van de maximum likelihoodschattingen behoorlijk ver-schillend. Dit is te zien in figuur 8. De theoretische verdeling is geen goede beschrijving van de empirische verdeling. In de tweede plaats geeft de methode van de probability weighted moments een schatting voor ξ die betekent dat deze methode per definitie niet betrouwbaar is. Hieruit volgt dat een rechtstreekse schatting van GEV-parameters door middel van de maximum-likelihood of via de methode van de probability weighted moments geen goede resultaten opleveren voor de reeks van dagelijkse neerslaghoeveelheden.
Nu volgen de resultaten van de pogingen een POT-model te schatten. De grafiek van de mean-excessfunctie geeft een stijgende rechte lijn te zien, met aan het begin een knik. Dit wordt veroorzaakt door de grote hoeveelheden waarnemingen van (net boven) 0. Een dui-delijke keuze voor u geeft de grafiek van de mean-excess functie zoals verwacht niet. Op grond van figuur 9 kan de grens overal neergelegd worden boven de 10.
Methode ξˆ βˆ σˆξ σˆβ xˆ0,95
MOM 0,10 85,2 0,2051 25,1 647
PWM 0,03 79,5 0,24 25,9 599
PKD -0,14 78,2 NaN NaN 541
Pickands-schatter geeft parameterschattingen die de data het best lijken te beschrijven, de momentenschatter geeft echter ook een goede ’fit’. De laatste geeft bovendien een schatting voor de standaarddeviatie. De standaarddeviatie is dusdanig dat moeilijk vastgesteld kan worden of de parameter ξ ongelijk is aan nul.
Zoals gesteld heb ik ook model geschat wat niet uitgaat van een identieke verdeling van de neerslaghoeveelheid voor iedere dag. Hier volgen de schattingsresultaten voor model-len met betrekking tot de jaarlijkse maxima van de neerslaghoeveelheid per etmaal. De resultaten van de recordanalyse zijn al gegeven. Verder blijkt niet uit de scatterplot dat er sprake is van afhankelijkheden. Een rechtstreekse maximum likelihoodschatting van de parameters van de GEV-verdeling geeft: ˆξ = −0, 18; ˆβ = 53, 1; ˆµ = 286, met bijbehorende standaarddeviaties van resp. 0,12; 7,0 en 8,2. De kwantielschatter, die een functie is van
deze parameters, is nu: ˆxp = 714. De kleine absolute waarde voor ξ geeft aanleiding om te
toetsen of hier sprake is van een gumbelverdeling. De zogenaamde gumbeltest geeft aan dat deze hypothese niet verworpen kan worden (overschrijdingskans is 35%). Een vergelijking tussen de empirische de theoretische verdeling geeft vertrouwen in dit model.
De mean-excess grafiek geeft min of meer een stijgende rechte lijn zien vanaf 250, zie figuur 11 De schattingsresultaten voor de parameters blijken ook redelijk constant vanaf deze grens. Een grens van 250 mm neerslag is echter moeilijk extreem te noemen, 88% van de waarnemingen zijn groter dan deze grens. De ’tussentijden’ zijn dan ook niet exponentieel verdeeld. Leggen we de grens wederom bij 350 dan blijkt dit laatste wel het geval te zijn. De schattingen voor het GPD-model zijn nu:
Methode ξˆ βˆ σˆξ σˆβ
MOM 0,10 85,2 0,2051 25,1
PWM 0,03 79,5 0,24 25,9
PKD -0,14 78,2 NaN NaN
Temperatuur
In deze paragraaf staan de resultaten van de pogingen om de gemiddelde dagelijkse tem-peraturen te modelleren. De temtem-peraturen zijn genoteerd in 0,1 graden celsius. In tegen-stelling tot de situatie bij neerslaggegevens is het bij temperatuurgegevens ook interessant extremen te modelleren die heel klein zijn. Hier geldt namelijk geen natuurlijke ondergrens van nul. Ik bespreek eerst de resultaten voor hoge extremen.
veel te hoog. Een simpele manier om deze moeilijkheid te omzeilen is het jaarlijkse maxi-mum of gemiddelde te berekenen en deze waarden te modelleren, evenals bij de neerslag. Een andere manier om het probleem van niet-identiek verdeelde variabelen op te lossen is het modelleren van de gecorrigeerde dagwaarden. De gemiddelde temperatuur in juni is anders dan in december, door het specifieke maandgemiddelde van iedere waarneming af trekken kan ik de waarnemingen normaliseren. Beide methodes heb ik onderzocht.
Het berekenen van de jaarlijkse maxima resulteert in 56 waarnemingen, zie figuur 12. Iedere waarneming beschrijft de maximale gemiddelde temperatuur die in het correspon-derende jaar gevonden is. De recordanalyse geeft geen blijk van niet-identiek verdeelde data: 5 records tegen een verwachting van 4,6 en bij een variantie 3,0. Een maximum
li-kelihoodschatting voor de parameters van de GEV-verdeling geeft: ˆξ = 0, 53; ˆβ = 20, 4 en
ˆ
µ = 234, 9. De kwantielschatter geeft vervolgens ˆx(0,95) = 37, 3. Op grond van figuur 13 stel
ik dat de Fr´echet-verdeling de waarnemingen goed beschrijft.
Het blijkt niet mogelijk op basis van deze waarnemingen een goed POT-model te bepalen. Het is moeilijk een grens te bepalen die bepaalt of een waarneming extreem is. De mean-excess plot geeft geen uitsluitsel en er is geen grens zodat de tijden tussen verschillende extreme waarnemingen exponentieel verdeeld zijn.
Nu volgen de resultaten van de analyse van de genormaliseerde reeks. Een histogram van de reeks genormaliseerde gemiddelde temperaturen staat in figuur 14. Een recordanalyse toegepast op deze reeks geeft geen bewijs tegen de veronderstelling dat deze bestaat uit identiek verdeelde variabelen. Er zijn 14 records waargenomen tegen een verwachting van 10.5 en een variantie van 4,8. Het schatten van de parameters van de GEV geeft:
Methode ξˆ βˆ µˆ σˆξ σˆβ σˆµ xˆ0.95
MLE 0,29 35,3 -11,8 0,02 0,18 0,24 157
PWM 0,28 34,1 -11.8 0,05 0,19 0,27 148
De ’fit’ blijkt heel goed te zijn, zie figuur 14. Deze reeks kan dus worden beschreven met een Fr´echet verdeling.
De mean-excess plot van deze reeks, zie figuur 15 geeft aan dat vanaf -5 de grafiek ongeveer rechtloopt. Vanaf een grens van 90 zijn de tussentijden tussen overschrijdingen bij benade-ring exponentieel verdeeld, getuige de qq-plot in figuur 16. De parameterschattingen voor de GPD hangen echter sterk af van de gekozen grens u. Het vertrouwen in dit model is daarom niet heel groot, ondanks het feit dat de GPD-verdeling de overschrijdingen boven de grens 90 redelijk benaderd.
bedragen ˆξ = 0, 28; ˆβ = 34, 1; ˆµ = −12, 05, de kwantielschatter ˆx0.95 = −21. De test
van gumbel geeft aan dat de vormparameter ξ daadwerkelijk groter is dan nul. Een POT model levert geen goede resultaten voor deze data. Dit wordt veroorzaakt doordat de mean-excessplot moeilijk is te interpreteren en de parameterschattingen sterk afhankelijk zijn van de gekozen grens u. Dit betekent praktisch gezien dat er teveel problemen zijn om een valide POT-model te maken voor de minimum waarnemingen.
Windsnelheid
De reeks die ik nu behandel bevat de maximale windsnelheid gemeten voor iedere dag tus-sen 1 januari 1951 en 31 december 2006. Een plot van de data staat in figuur 21, deze geeft een indicatie dat er sprake is van seizoenseffecten in de windsnelheid. Een recordanalyse geeft echter hier geen uitsluitsel over en een correctie voor het maandelijks gemiddelde maakt de uitslag van de toets niet duidelijker. De reeks met jaarlijkse maxima is zonder twijfel wel constant verdeeld over de tijd getuige ook de resultaten van de recordanalyse:
Reeks Verwacht Variantie Waargenomen
Dagdata max. windsn. 10,5 8,8 5
Gecorr. dag. max.windsn. 10,5 8,8 6
Jaarl. max. windsn. 4,6 3,0 6
Een scatterplot geeft geen sterke indicatie dat de data afhankelijk zouden zijn. Ik geef de resultaten van de pogingen een model te schatten voor de dagdata. Een rechtstreekse schatting van de parameters van de GEV levert een uitstekende ’fit’ van de data, zie fig 22. De waarden voor de parameterschattingen zijn:
Methode ξˆ βˆ µˆ xˆ0.95
MLE 0,11 40,4 87,4 190
PWM 0,09 38,3 87,1 192
De covariantiematrix bij deze schattingen is ongeveer constant en levert relatief kleine waarden, de standaarddeviaties komen uit op : ˆσξ = 0, 005, ˆσβ = 0, 21 en ˆσµ = 0, 30. Voor
de pwm-schatting volgen ongeveer dezelfde waarden. Op grond van de gumbeltest van gumbel, kan gesteld worden dat de parameters significant verschillen van 0. Ik merk hierbij op dat de overschrijdingskans dicht bij 5% zit. Dit feit gecombineerd met de vorm van een histogram van de data geeft mij het sterke vermoeden dat de data het beste gemodelleerd kunnen worden een Gumbel-verdeling.
De mean-excess functie geeft een rechte lijn vanaf ongeveer 100. De GPD-parameters zijn afhankelijk van de gekozen grens u, zie figuur 24. Vanaf 250 zijn de tussentijden exponen-tieel verdeeld, de parameterschattingen voor het GPD behorende bij deze grens staan in onderstaande tabel:
Methode ξˆ βˆ xˆ0.95
MOM 0,17 34,0 382
PWM 0,25 32,4 394
De standaarddeviaties verschillen niet heel sterk in absolute waarde. De kleinste standaard deviaties worden verkregen met de momentenschatter: ˆσβ = 0, 09 en ˆσµ= 4, 2. Dit betekent
op grond van de gumbeltest dat ξ niet significant verschilt van nul. Het verschil in afstand tussen de theoretische en de empirische verdeling lijkt het kleinst bij het resultaat van de momentenschatter. Een grafiek van deze twee verdelingen staat in figuur 25. De kwaliteit van het model lijkt heel acceptabel. Pickands’ schatter geeft geen goede beschrijving van de data, bovendien convergeren de schattingen voor de covariantiematrix niet. Pearsons chis-quare toets levert een overschrijdingskans van 1; een QQ-plot van de tussentijden tegen een willekeurige exponenti¨ele verdeling geeft alleen in de staarten aanleiding tot twijfel. Opvallend is dat de POT-modellen een hogere schatting opleveren voor het 0.95-kwantiel dan het model op basis van de GEV-verdeling. Omdat de schattingen voor het laaste model op meer waarnemingen zijn gebaseerd geef ik hieraan de voorkeur.
Het GEV-model voor de jaarlijkse maximum windsnelheid is wederom gebaseerd op 56 waarnemingen. Het schatten van een GEV-verdeling levert het volgende resultaat op:
Methode ξˆ βˆ µˆ xˆ0.95
MLE 0,11 40,4 87,4 342
PWM 0,05 32,5 257,1 347
De gumbeltest geeft aan dat de parameter ξ niet significant van nul verschilt (overschrij-dingskans 33%), ik concludeer dat we te maken hebben met een Gumbel-verdeling. De kwaliteit van het model is op basis van figuur 23 goed te noemen.
Uren Zonneschijn
Nu volgen de resultaten van het toepassen van extreme waarden theorie op de reeks over uren zonneschijn per etmaal. Deze reeks bestaat uit relatief veel waarnemingen gelijk aan 0, corresponderend met de dagen waarop de zon niet schijnt. Dit is de natuurlijke ondergrens voor de waarden in deze reeks. Het aantal uren zon blijkt na aanleiding van de recordana-lyse niet identiek verdeeld over het jaar. Dit is logisch, in de winter komt de zon later op en gaat deze eerder onder. Na verwijdering van de seizoensinvloeden zijn er nog steeds problemen met de stabiliteit van de verdeling. Bovendien geeft de scatterplot aan dat er sterke afhankelijkheden verstopt zitten in de data. Waarschijnlijk heeft dit te maken met het probleem dat ook binnen een maand nog verschillen zitten in de opkomst en ondergang van de zon, maar bewijs heb ik hier niet voor gevonden.
Reeks Verwacht Variantie Waargenomen
Ongecorrigeerde Dagdata 10,5 8,8 20
Gecorrigeerde Dagdata 10,5 8,8 20
Jaargemiddelde 4,6 3,0 4
Een toets op onafhankelijkheid d.m.v de scatterplot van de serie geeft geen aanleiding tot twijfel, daarom probeer ik van deze reeks de extremen te modelleren. Een plot van de data leert dat er een paar extremen te zien zijn, maar het aantal is niet heel groot. Ik probeer eerst of deze waarnemingen kunnen worden gemodelleerd met de GEV. Parameterschat-tingen met behulp van de maximum likelihood (ML) en de probability weighted moments (PWM) geven het volgende resultaat:
Methode ξˆ βˆ µˆ xˆ0.95
MLE 0,09 3,93 40.1 50.3
PWM 0,11 4,6 40.1 52,0
De ’fit’ is redelijk te noemen, alleen in de linker staart lijkt de geschatte verdeling wat min-der goed overeen te komen dan met de empirische verdeling, zie figuur 18. De mean-excess functie behorende bij deze reeks staat in figuur 20. Deze loopt recht vanaf ongeveer 36 . Vanaf de grens 40 kan gesteld worden dat de tussentijden exponentieel verdeeld zijn. Hierbij moet opgemerkt worden dat dit niet heel duidelijk vastgesteld kan worden, de qq-plot laat ruimte voor twijfel vanwege afwijkingen in de extreme staarten. De schattingsresultaten voor de parameters van de GPD-verdeling vermeld ik hieronder:
Methode ξˆ βˆ βˆ µˆ xˆ0.95
MOM 0,32 6,02 0,18 1,39 66
PWM 0,44 6,56 0,24 1,66 77
PKD 0,3884 5,4057 NaN NaN 66
De momentenschatter heeft de kleinste standaarddeviatie en de beste ’fit’ van de data, zie figuur 19.
Luchtdruk
Reeks Verwacht Variantie Waargenomen
Ongecorrigeerde dagluchtdruk 10,5 8,8 16
Gecorrigeerde dagluchtdruk 10,5 8,8 20
Gemiddelde maandelijkse luchtdruk 7,1 5,4 6
Allereerst volgen de resultaten van het modelleren van de maxima. Het toepassen van de GEV-verdeling op de reeks met maandelijkse gemiddelden geeft geen bevredigende resulta-ten. Dit wordt veroorzaakt door de bijzondere vorm van de empirische verdeling, zie figuur
276. Daarom probeer ik verder te gaan met een POT-model.
De mean-excessfunctie is ongeveer recht vanaf 1045 millibar, de tussentijden zijn echter nooit exponentieel verdeeld. De schatters blijken wederom een functie van de grens u (fi-guur 28). Het schatten van het POT-model voor een grens van 1045 millibar levert desalnie-temin een goede theoretische verdeling op ten opzichte van de empirische verdelingsfunctie van de GPD. De resultaten van de schatters staan in de tabel:
Methode ξˆ βˆ βˆ µˆ xˆ0.95
MOM 0,49 0.90 0.12 0.11 1700
PWM 0,59 0.96 0,15 0.13 1835
PKD 0,66 0.88 NaN NaN 1874
De schatters zijn echter sterk afhankelijk van de grens u, dit geeft geen vertrouwen in het model.
Voor het modelleren van minima van de maandelijkse reeks gemiddelde luchtdruk gelden eigenlijk dezelfde resultaten. De GEV-verdeling levert geen goede beschrijving op van de reeks. Het schatten van een POT-model is nauwelijks succesvoller. De meanexcess-functie is weliswaar recht vanaf 9,85, maar de tussentijden zijn echter niet exponentieel verdeeld. Daarnaast zijn ook hier de parameterschattingen sterk afhankelijk van de grens u.
Resultaten familieselectie en parameterschattingen
In de vorige paragraaf heb ik de resultaten gepresenteerd van de analyses van de afzonder-lijke reeksen. In deze paragraaf presenteer ik de resultaten die ik gevonden naar aanleiding van een analyse van de samenhang tussen verschillende reeksen waarnemingen. De basis van deze analyse vormen de modellen die ik gemaakt heb in de vorige paragraaf. Voor iedere reeks heb ik geprobeerd marginale verdeling te schatten, al dan niet na de nodige transformaties. Om vervolgens een verband te kunnen modelleren tussen twee verschillen-de marginale ververschillen-delingen heb ik verschillen-de betreffenverschillen-de waarnemingen gekoppeld. Ieverschillen-dere koppeling is gebasseerd op het tijdstip waarop de waarnemingen gedaan zijn. Bijv. de gemiddelde
temperatuur op 1 januari 1980 heb ik gekoppeld aan de neerslagsom op 1 januari 1980, de gemiddelde temperatuur op 2 januari 1980 aan de neerslagsom op 2 januari 1980 enz. Dit is een belangrijke reden waarom ik in het vorige hoofdstuk gebruik gemaakt heb van reeksen van gelijke lengte. Een gevolg van deze koppeling is dat niet alle modellen verge-leken kunnen worden. Het model voor het jaarlijks maximum van de hoeveelheid neerslag
die op ´e´en dag gevallen is kan niet vergeleken worden met het model voor de dagelijkse
hoeveelheid neerslag.
Ik heb gebruik gemaakt van de methode van Genest en de variant gebaseerd op de
mini-malisatie van de L2-norm. De methode van Genest en Rivest heb ik alleen gebruikt om te
selecteren tussen de Gumbel- en de Frankcopulafamilie. De andere families gaven te gecom-pliceerde parameterschattingen. Voor de variant op Genest en Rivest heb ik een selectie gemaakt tussen de familietypen Clayton, Gumbel, Frank en Log-copula. Er bestaan meer familietypen, maar om het onderzoek niet te groot te laten worden heb ik me beperkt tot deze typen.
Ik heb alle modellen uit het vorige hoofdstuk op samenhang onderzocht. In dit hoofdstuk presenteer ik de resultaten van de modellen die tot een succesvolle beschrijving van de sa-menhang in reeksen leiden. Dit laatste meet ik met de toetsgrootheid ’Pearsons Goodness-Of-Fit’.
Het verband tussen zonneschijn en temperatuur
Het gebruik van copula’s kan uitstekend gedemonstreerd worden aan de hand van een model die de samenhang tussen zonneschijn en temperatuur beschrijft. Dat er een verband bestaat tussen beide grootheden kan worden beschouwd als ’Commen sense’. Daarom leent deze situatie zich uitstekend voor een demonstratie van de behandelde methoden. In het vorige hoofdstuk heb ik beide reeksen individueel gemodelleerd. Hieruit volgde dat het jaarlijkse temperatuurgemiddelde gemodelleerd kan worden met de GEV-verdeling met parameters ξ = 0, 38 ; β = 66, 3 en µ = 77, 8 en dat het jaarlijkse gemiddelde aantal uren zon gemodelleerd kan worden door een GEV-verdeling met parameters ξ = 0, 11; β = 4, 6 en µ = 40, 1.
Ik koppel nu de gemiddelde temperatuur in een jaar aan het gemiddeld aantal uren zon in datzelfde jaar. Ik heb de volgende 56 waarnemingen waarnemingen als resultaat:
(X1, Y1), (X1, Y1), . . . , (X56, Y56)
Iedere component heeft bovendien een bekende marginale verdeling :
Xi ∼ GEV (0, 38; 66, 3; 77, 8), i = 1, . . . , 56.
Yi ∼ GEV (0, 11; 4, 6; 40, 1), i = 1, . . . , 56.
Ik ben nu op zoek naar een copulafunctie die de bivariate verdeling FX,Y(x, y) weergeeft.
Figuur 1: Grafische weergave procedure Genest en Rivest voor zon en temperatuur.
lijn in het midden is het empirische Kendalls proces Kn(v). Deze is geconstrueerd op
ba-sis van vergelijking 5. De twee lijnen die hieraan min of meer evenwijdig lopen vormen de betrouwbaarheidsintervallen, deze zijn gevonden met behulp van vergelijking 6. Voor de Claytonfamilie zijn de betrouwbaarheidsintervallen exact, voor de Gumbelfamilie een benadering.
De gele continue lijn is het theoretische ’Kendalls process’ op basis van de Claytons familie met parameter α = 0, 65. Deze parameterwaarde is bepaald op basis van de vergelijking
ˆ
α = 2−ˆτˆτ. De empirische waarde van Kendalls tau, ˆτ die iets zegt over de mate van sa-menhang tussen beide reeksen heeft een waarde van 0,39. Dit betekent dat er een positief verband bestaat tussen beide reeksen, hetgeen overeenkomt met hetgeen verwacht kan wor-den. De gestreepte lijn is het theoretische ’Kendalls process’ op basis van Gumbels familie met parameter α = 1, 3. Deze parameterwaarde is bepaald op basis van de vergelijking
ˆ
α = 1−ˆτˆτ. Uit de figuur blijkt dat de theoretische waarden niet altijd binnen het
95%-betrouwbaarheidsinterval vallen en lijken daarom tekort te schieten in de beschrijving van het empirische proces.
Op grond van dit plaatje lijkt de copulafunctie uit de Claytons familie met α = 0, 65 nog het beste in staat de data te beschrijven. Ik zal het model nu formeel toetsen met Pearson’s Goodness-Of-Fit test. De matrix met geobserveerde waarden is:
Y19 Y37 Y56
X19 10 6 2
X37 6 8 4
X56 2 4 12
In de middelste cel staat dus het aantal waarnemingen (Xi, Yi) waarvoor geldt: X18 < Xi <
de volgende waarden verwacht worden: y0,33 y0,66 y1
x0,33 9, 9 5, 2 3, 6
x0.66 5, 2 6, 8 6, 7
x1 3, 6 6.7 8, 4
In bovenstaande matrix staat y.33 voor het 33%-kwantiel van de marginale verdeling van
gemiddelde zonneschijn per dag. De GOF-toets geeft een overschrijdingskans van 70%. Een sterke indicatie dat het model een goede beschrijving biedt van de waarnemingen. De GOF-toets voor het model op basis van Franks familie met parameter α = 1, 3 geeft de volgende matrix met verwachte waarden:
y0,33 y0,66 y1
x0,33 12, 7 5, 1 0, 9
x0.66 5, 1 9, 4 4, 1
x1 0, 9 4, 1 13, 7
De bijbehorende overschrijdingskans is 81%. Hierbij heb ik cellen met minder dan 5 waarne-mingen bijelkaar genomen. Beide copulafuncties lijken geschikt om de data te beschrijven. De variant op de procedure van Genest en Rivest levert een nog overtuigender model. Zo-als in figuur 2 te zien is, beschrijven de copulafuncties de data nog veel beter. Dit is geen verrassend resultaat, want de ruimte voor het kiezen van de juiste functie is veel groter omdat uit meer parameters gekozen kan worden. De rode lijn geeft de beste theoretische
copulafunctie weer. Dat is de copulafunctie die op grond van de discrete L2 het dichtste
bij de waargenomen gele lijn ligt, in dit geval de log-copula met parameters α = 0,1672 en γ = 0,1507. De gele ’schokkerige’ lijn beschrijft het empirische ’Kendall’s process’, de twee evenwijdige groene lijnen de 95%-betrouwbaarheidsintervallen. Verder geeft de en-kelvoudig gestreepte blauwe lijn de beste copulafunctie uit de Claytonfamilie, de blauwe gestippelde lijn de beste copulafunctie uit de Frankfamilie en de afwisselend gestreepte en gestippelde blauwe lijn de beste copulafunctie uit de Gumbelfamilie. De log-copulafunctie CL(0, 17; 0, 15) leidt tot de volgende verwachte waarden:
y.33 y.66 y1
x.33 11, 1 5, 8 1, 8
x.66 5, 8 8, 1 4, 8
x1 1, 8 4, 8 12, 1
De GOF-toets geeft hier een overschrijdingskans van 99,99%.
Figuur 2: Grafische weergave procedure Genest en Rivest
Neerslag en luchtdruk
Voor de maandelijkse gemiddelde neerslag per dag heb ik in een GEV-verdeling gevonden
met parameters ˆξ = −0, 18; ˆβ = 53, 1; ˆµ = 286. Voor maandelijkse gemiddelde luchtdruk
vond ik geen betrouwbaar model. Toch kan er een goede samenhang gevonden worden tus-sen deze reektus-sen. Er doet zich nog een opvallend feit voor. De variant op de procedure van Genest en Rivest selecteert een copulafunctie die het empirische Kendalls proces uitstekend beschrijft getuige figuur 5. De gele ’schokkerige’ lijn beschrijft het empirische ’Kendall’s process’, de twee evenwijdige groene lijnen de 95%-betrouwbaarheidsintervallen. Verder geeft de enkelvoudig gestreepte blauwe lijn de beste copulafunctie uit de Claytonfamilie, de blauwe gestippelde lijn de beste copulafunctie uit de Frankfamilie en de afwisselend gestreepte en gestippelde blauwe lijn de beste copulafunctie uit de Gumbelfamilie. Een GOF-toets met 9 cellen geeft een overschrijdingskans van 95%. Omdat het aantal waarne-mingen per cel redelijk groot is, breid ik het aantal cellen uit naar 16. Op dat niveau zien geeft de GOF-toets een overschrijdingskans van minder dan 1%. De onderstaande matrices geven de waargenomen en verwachte waarden met 9 cellen:
