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.
De eerste methode stoelt op de veronderstelling dat de data afkomstig zijn (uit het attrac-tiedomein van) de GEV-verdeling. Op basis van deze aanname zal ik gebruik maken van schatters die de parameters van de GEV-verdeling bepalen. Ik zal vooral gebruik maken van de maximum likelihoodschatter en de methode van de probability weighted moments.
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.
Hoewel de recordanalyse geen redenen tot twijfel geeft, ben ik op praktische gronden niet geheel overtuigd van de juistheid van laatstgenoemde veronderstelling. Het is goed denk-baar dat een grote neerslaghoeveelheid op dag 1 een gevolg heeft voor de neerslagkans op de volgende dag. Bovendien kunnen seizoensinvloeden invloed hebben op de neerslagkans. Om mogelijke afhankelijkheden en seizoensinvloeden in de data uit te sluiten heb ik beslo-ten ook alle jaarlijkse maxima van de neerslaghoeveelheid per etmaal modelleren. Deze zijn in ieder geval identiek verdeeld, klimaatveranderingen daargelaten. Voor deze reeks heb ik ook een recordanalyse gedaan, een samenvatting van deze en de vorige recordanalyse staat
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.
De keuze voor u is echter wel een belangrijke, de schattingen voor de parameters ξ en β zijn niet onafhankelijk van de keuze voor u, zie grafiek 10. De schattingen, gebaseerd op de benadering door de paretoverdeling, voor ξ liggen rond de nul. Voor een lage grens liggen ze er meestal onder, voor een hoge grens erboven, zie de grafiek 10. Voor β lopen ze uiteen van 20 to 80, van 20 tot 60 min of meer evenredig met de grenswaarde u. De variantie van de schatters loopt vanzelfsprekend op naarmate de grens hoger wordt aangezien er minder waarnemingen gebruikt worden voor de schatting. Een andere leidraad bij de keuze van de grens u is de verdeling van de overschrijdingen. Wanneer deze boven een bepaalde grens Poisson verdeeld zijn, is dat een indicatie voor de juiste grens. Ik toets dit aan de hand van de tussentijden. Als de overschrijdingen arriveren volgens een Poisson process zullen deze exponentieel verdeeld zijn. De qq-plot geeft voor de overschrijdingen over 350 een goede indicatie dat de tussentijden inderdaad exponentieel verdeeld zijn en overschrijdin-gen voloverschrijdin-gens een Poisson proces arriveren. Bij deze grens geeft het POT-model de voloverschrijdin-gende schattingsresultaten voor de parameters van de GPD-verdeling en een kwantiel-schatter:
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.
Voorafgaand aan de analyse is het het goed om te constateren dat op een dag in juni een andere temperatuur verwacht wordt dan op een dag in december. Het zou onzin zijn om te veronderstellen dat de waarnemingen identiek verdeeld zijn. Het aantal records blijkt ook
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.
Nu volgen de resultaten van het modelleren van de lage extremen, de minima. Voor het modelleren van de minima heb ik de data op dezelfde manier gecorrigeerd als bij de maxi-ma: Van iedere waarneming heb ik het betreffende maandgemiddelde afgetrokken. Er zijn 15 records waargenomen tegen een verwachting van 10.5 en een variantie van 4,8. Een probability weighted moments schatter geeft een goede ’fit’ met de data (figuur 17), dit in tegenstelling tot de maximum likelihood schatter. De schattingen voor de parameters
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:
5. De negatieve locatie parameter wordt veroorzaakt doordat ik het negatief van de waarnemingen geno-men heb.
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.
Op grond van deze redenen heb er ik ervoor gekozen alleen een model te maken voor het jaargemiddelde, de totale hoeveelheid zon in een jaar gedeeld door het aantal dagen. Bij deze datareeks van 56 waarnemingen blijkt geen bewijs tegen de hypothese dat de jaarge-middelden gelijk en onafhankelijk verdeeld zijn.
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