FACULTEIT WETENSCHAPPEN Vakgroep Toegepaste Wiskunde, Informatica en Statistiek Temperatuurderivaten Sarah Vanhuylenbroeck Promotor: Prof. dr. M. Vanmaele Masterproef ingediend tot het behalen van de academische graad van master in de wiskunde, afstudeerrichting toegepaste wiskunde. Academiejaar 2013 - 2014 Voorwoord Door te kiezen voor de minor economie in mijn bachelor en de minor economie en verzekeringen in mijn masteropleiding, wilde ik in mijn masterproef een onderwerp behandelen dat bij deze opleidingen aanleunt. De keuze voor een onderwerp in verband met financiele wiskunde was dus snel gemaakt. Ik wilde ook graag een onderwerp behandelen dat aansluit bij het dagelijkse leven. Toen mijn promotor Prof. dr. M. Vanmaele mij vertelde dat recent de focus eerder op weerderivaten, zoals afgeleide producten met wind en temperatuur als onderliggende, gevestigd is en daar de aandacht in het bijzonder gaat naar het modelleren van wind en temperatuur, heb ik besloten om mijn masterproef aan temperatuurderivaten te wijden. Het schrijven van deze masterproef was een zeer leerrijke ervaring en ik ben dan ook erg blij dat ik voor dit onderwerp gekozen heb. Eerst en vooral wil ik mijn promotor Prof. dr. M. Vanmaele bedanken voor de tijd die zij vrijgemaakt heeft om mij te begeleiden bij het schrijven van deze masterproef. Ik wil haar ook danken voor de raad die ze mij gegeven heeft en voor het nalezen en verbeteren van eerdere versies van mijn masterproef. Ik dank ook het KMI, en in het bijzonder Alex Dewalque, omdat zij hun temperatuurgegevens ter beschikking stelden voor mijn masterproef. Vervolgens wil ik ook mijn ouders bedanken, omdat zij het voor mij mogelijk hebben gemaakt om te studeren en zij nog steeds op alle mogelijke manieren voor mij klaarstaan. Verder bedank ik ook Linda Golden, Chuanhou Yang, Maximilian Wimmer en David Stillberger omdat zij de tijd namen om mijn vragen over hun papers te beantwoorden. Tot slot dank ik ook mijn vrienden en mijn vriend, omdat zij voor de nodige ontspanning zorgden, niet alleen tijdens de periode waarin ik mijn masterproef schreef, maar in mijn volledige carri`ere als student. i ii Toelating tot bruikleen De auteur geeft de toelating deze masterproef voor consultatie beschikbaar te stellen en delen van de masterproef te kopi¨eren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze masterproef. Sarah Vanhuylenbroeck, mei 2014 iii iv Inhoudsopgave Inleiding 1 1 Temperatuurderivaten en temperatuurindices 3 2 Modellering van de temperatuur 7 2.1 Dynamische econometrische modellen . . . . . . . . . . . . . . . . . . . . 7 2.1.1 MA-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.2 AR-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.3 ARMA-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.4 ARIMA-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.5 SARIMA-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.6 GARCH-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 CAR-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Tijdscontinue modellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 Ornstein-Uhlenbeck . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Hull-White . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Praktijk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.1 Het model van Campbell en Diebold . . . . . . . . . . . . . . . . 19 2.4.2 Het Alaton-model . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4.3 Het Benth-model . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4.4 Econometrische modellen voor Belgische temperatuurdata . . . . 29 Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.4 2.5 3 Temperatuurderivaten prijzen 55 3.1 Black-Scholes-Merton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.2 Historical burn analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.2.1 Prijzen op basis van HBA met Belgische temperatuurdata . . . . 62 Indexmodellering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 v 3.4 3.5 3.6 Dagelijkse simulatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.4.1 Discreet proces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.4.2 Continu proces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.4.3 Prijzen op basis van discrete dagelijkse simulatie . . . . . . . . . . 70 3.4.4 Prijzen op basis van continue dagelijkse simulatie . . . . . . . . . 71 3.4.5 Prijzen met behulp van weersvoorspellingen . . . . . . . . . . . . 97 Alternatieve methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.5.1 Indifference valuation . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.5.2 Equilibrium valuation . . . . . . . . . . . . . . . . . . . . . . . . 111 Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 A Anderstalige samenvatting 119 B R-code 123 B.1 Algemene eigenschappen van de tijdreeks . . . . . . . . . . . . . . . . . . 123 B.2 Differentie-operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 B.3 Fourierreeks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 B.4 Prijzen op basis van discrete dagelijkse simulatie . . . . . . . . . . . . . . 133 C Figuren 135 vi Inleiding Op klimatologisch vlak leek 2013 een jaar van grote extremen. De koude winter, met nog een sneeuwstorm op 12 maart, werd moeiteloos gevolgd door een prachtige zomer. Dit soort onverwachte weersomstandigheden houden echter een zeker risico in voor bepaalde bedrijven (denk bijvoorbeeld aan bioscopen die tijdens een mooie zomer het aantal bezoekers zien dalen) en hebben hierdoor een grote invloed op de economie. Dit betekent dat het weer steeds vaker een bijkomende zorg wordt voor de risicomanagers van bedrijven. De financi¨ele markt biedt sinds 1996 − 1997 echter een product aan om weerrisico’s te hedgen, namelijk weerderivaten. Deze speciale vorm van derivaten is er vooral gekomen onder impuls van de elektriciteitsbedrijven. Zij zijn immers zeer gevoelig voor de weersomstandigheden (gedurende een zachte winter verwarmen de mensen hun huizen minder) en zochten naar een middel om deze gevoeligheid te milderen. Weerderivaten zijn geen verzekeringen, dit wil zeggen dat ze ingevoerd werden om niet-catastrofale weersomstandigheden te dekken. Op die manier zullen zij niet de schade door een orkaan (dit komt zelden voor) dekken, maar wel bijvoorbeeld het verlies aan inkomsten van een horecazaak aan de kust door een slechte zomer. Derivaten, in het algemeen, zijn financi¨ele instrumenten die hun waarde ontlenen aan een onderliggend actief (de onderliggende), bijvoorbeeld aandelen of goud. Weerderivaten hebben als specifieke onderliggende het weer (met name temperatuur, regen, sneeuwval,. . . ). Zoals de andere financi¨ele derivaten kunnen weerderivaten zowel gebruikt worden om een risico te hedgen als voor speculatie. Er is echter ook een belangrijk verschil met andere derivaten. De onderliggende (d.i. het weer) is niet verhandelbaar en heeft hierdoor geen waarde. Om dit probleem op te lossen worden weerindices ingevoerd, deze worden dan gebruikt als onderliggend actief. De meerderheid van de verhandelde weerderivaten wordt geschreven op temperatuurindices. Dit is hoofdzakelijk het geval omdat de meerderheid van de bedrijven die in weerderivaten handelen energiebedrijven zijn en hun inkomsten vooral afhankelijk zijn van de temperatuur. Hierdoor zullen we ons in deze studie dan ook beperken tot temperatuurderivaten. Deze studie is tweeledig, zowel het modelleren van 1 2 de temperatuur als het prijzen van een derivaat zal worden besproken. Deze masterproef is opgebouwd uit drie hoofdstukken. Hoofdstuk 1 is eerder een inleidend hoofdstuk waarin we de definities geven van de temperatuur en de verschillende temperatuurindices die we beschouwen. Deze definities zullen doorheen de masterproef gebruikt worden. In Hoofdstuk 2 behandelen we het eerste aspect van deze studie, namelijk het modelleren van de temperatuur. Eerst en vooral belichten we de theorie van de vaak gebruikte econometrische basismodellen en stochastische differentiaalvergelijkingen. Daarna beschrijven we in Sectie 2.4 enkele modellen die in de literatuur vaak besproken worden. We stellen, tot slot, ook zelf twee econometrische modellen voor om de Belgische temperatuur te beschrijven. Het eerste model specifieert de seizoenseffecten aan de hand van de differentie-operator, terwijl het tweede model een Fourierreeks gebruikt om deze effecten voor te stellen. Dit laatste model gaan we vervolgens in Hoofdstuk 3 gebruiken om een fictieve calloptie te prijzen. Het derde en laatste hoofdstuk behandelt het tweede aspect van deze studie, het effectief prijzen van temperatuurderivaten. We bespreken er verschillende methoden om temperatuurderivaten te gaan prijzen. Eerst en vooral gaan we kijken naar het BlackScholes-Merton-model, maar dit bekende model voor het prijzen van financi¨ele derivaten kan voor temperatuurderivaten niet gehanteerd worden. Hierdoor bestaan er heel wat alternatieven. De eenvoudigste van deze methoden is historical burn analysis en deze aanpak wordt dan ook als eerste besproken. Aangezien deze methode simpel uit te voeren is, berekenen we in Sectie 3.2.1 de numerieke prijs van een fictieve calloptie gebruik makend van Belgische temperatuurdata. Een volgende, iets complexere methode is indexmodellering. Deze methode onderstelt een bepaalde verdeling voor de gebruikte temperatuurindex. Daaropvolgend bespreken we een techniek die nauwkeuriger is dan de twee voorgaande methodes. Dit is het prijzen van temperatuurderivaten aan de hand van dagelijkse simulatie. Deze methode maakt gebruik van een model voor de temperatuur om een derivaat te prijzen, hier komt de kennis uit Hoofdstuk 2 dus van pas. We bespreken in Sectie 3.4.3 en 3.4.4, respectievelijk, hoe we een discreet en een continu model voor de temperatuur kunnen gebruiken om temperatuurderivaten expliciet te gaan prijzen. Tot slot bespreken we in Sectie 3.5 enkele methoden die ontwikkeld werden om derivaten in incomplete markten te prijzen. Omdat de markt voor het weer incompleet is, kunnen deze methoden ook gebruikt worden om temperatuurderivaten te prijzen. Het zijn voornamelijk methoden die gebruik maken van zogenaamde nutsfuncties. Hoofdstuk 1 Temperatuurderivaten en temperatuurindices Zoals reeds in de inleiding vermeld werd, zijn temperatuurderivaten financi¨ele instrumenten waarvan de payoff afhankelijk is van de temperatuur. Daar temperatuur niet kan worden opgeslagen en dus ook niet verhandelbaar is, heeft deze geen prijs. Omdat dit betekent dat de onderliggende geen waarde heeft, worden vaak zogenaamde temperatuurindices gebruikt als onderliggende. Volgens Alexandridis and Zapranis (2013) zijn de meest gebruikte indices de Heating Degree Days (HDD), de Cooling Degree Days (CDD), de Cumulative Average Temperature (CAT) en de Pacific Rim (PR) index. Vooraleer we deze indices kort bespreken, merken we nog op dat de temperatuur volgens Alaton et al. (2002) als volgt wordt gedefinieerd. Definitie 1. Gegeven een specifiek weerstation. Zij Timax en Timin respectievelijk de maximum- en minimumtemperatuur (in graden Celsius) gemeten op dag i, dan is de temperatuur op dag i Ti = Timax + Timin . 2 (1.1) De meerderheid van de temperatuurderivaten is gebaseerd op de accumulatie van HDDs of CDDs tijdens een bepaalde periode (een maand of de winter/zomerperiode) en deze worden het vaakst gebruikt in de VSA, Canada en Australi¨e. Gewoonlijk omvat het HDD-seizoen de wintermaanden van november tot en met maart en omvat het CDDseizoen de zomermaanden van mei tot en met september; april en oktober worden vaak de flankmaanden genoemd. De referentiewaarde voor de HDD- en CDD-index is 18◦ C. De reden hiervoor is dat 3 4 Hoofdstuk 1. Temperatuurderivaten en temperatuurindices mensen bij een temperatuur lager dan 18◦ C hun huis verwarmen (heating) en bij een temperatuur boven 18◦ C eerder hun airco gebruiken (cooling). Een HDD staat dan ook voor het aantal graden dat de dagtemperatuur lager is dan 18◦ C en een CDD geeft het aantal graden aan dat de dagtemperatuur boven 18◦ C is. Samengevat krijgen we de volgende definitie volgens Alaton et al. (2002). Definitie 2. Zij Ti de temperatuur op dag i. We defini¨eren de heating degree days, HDDi , en de cooling degree days, CDDi , gegenereerd op die dag als HDDi = max{18 − Ti , 0}, (1.2) CDDi = max{Ti − 18, 0}, (1.3) en respectievelijk. In Europa worden hoofdzakelijk temperatuurderivaten voor de zomer met als onderliggende de CAT-index verhandeld. De CAT-index wordt berekend als de som van de gemiddelde dagtemperaturen tijdens de contractperiode. Zoals in Benth and Benth (2012), nemen we als referentie voor de gemiddelde dagtemperatuur het gemiddelde van de maximum- en minimumtemperatuur, zoals in Definitie 1. We krijgen dan de volgende definitie voor de CAT-index. Definitie 3. Zij Ti de temperatuur op dag i, dan wordt de cumulative average temperature index over een periode [t1 , t2 ] gedefinieerd als CAT(t1 , t2 ) = t2 X Ti , (1.4) i=t1 waar de periode meestal een kalendermaand of een bepaald seizoen is. De Pacific Rim index is nog een ander soort temperatuurindex, die hoofdzakelijk gebruikt wordt in Azi¨e, met name Japan. In Alexandridis and Zapranis (2013) wordt deze index simpelweg gedefinieerd als het gemiddelde van de CAT-index over een specifieke periode. Definitie 4. Zij CAT(t1 , t2 ) de CAT-index bepaald over de periode [t1 , t2 ], dan wordt de Pacific Rim index als volgt gedefinieerd Pt2 PR(t1 , t2 ) = Ti CAT(t1 , t2 ) = . t2 − t1 t2 − t1 i=t1 (1.5) 5 De gebruikte temperatuurindex is een belangrijke factor, het is namelijk ´e´en van de parameters die het contract definieert. Volgens Schiller et al. (2012), Alaton et al. (2002) en Alexandridis and Zapranis (2013) worden temperatuurderivaten bepaald door de volgende parameters: • de onderliggende index; • het type contract (calloptie, putoptie, future,. . . ); • de looptijd (contractperiode); • het uitoefenniveau; • de tick size, d.i. de hoeveelheid geld die de houder van het contract ontvangt voor elke indexwaarde boven het uitoefenniveau, vaak gedefinieerd als een bedrag per indexpunt; • het weerstation waarvan de gegevens worden verkregen; • de maximale uitbetaling (als er een bovengrens voor uitbetaling is); • de premie die door de koper aan de verkoper betaald wordt (bespreekbaar). 6 Hoofdstuk 1. Temperatuurderivaten en temperatuurindices Hoofdstuk 2 Modellering van de temperatuur Omdat temperatuurderivaten de temperatuur(index) als onderliggende hebben, zullen we in dit hoofdstuk een model proberen te vinden dat de temperatuur beschrijft. In de literatuur worden doorgaans twee groepen van modellering beschouwd, zoals in Schiller et al. (2012) onderscheiden we de aanpak die gebruik maakt van econometrische (discrete) processen enerzijds en de tijdscontinue aanpak anderzijds. Beiden hebben het voorspellen van de temperatuur a.d.h.v. een geconstrueerd proces, dat door het model beschreven wordt, tot doel. De econometrische aanpak is gesteund op het bouwen van discrete, econometrische modellen die de eigenschappen van de temperatuur implementeren. De tijdscontinue aanpak steunt op stochastische differentiaalvergelijkingen die de temperatuursveranderingen beschrijven. Omdat de uiteindelijke prijs van het derivaat deels bepaald wordt door de temperatuur(index), is het belangrijk dat het model de realiteit goed weergeeft. Zo zal het verkregen model rekening moeten houden met de seizoenen, aangezien de temperatuur in de winter wel degelijk lager is dan in de zomer. Er bestaan verschillende modellen die de temperatuurschommelingen beschrijven, velen zijn echter slechts een toepassing/aanpassing van onderstaande basismodellen. 2.1 Dynamische econometrische modellen In de literatuur zijn er verschillende studies terug te vinden waarin de temperatuur gemodelleerd wordt aan de hand van econometrische modellen. Een mooi voorbeeld hiervan is het model dat voorgesteld wordt door Campbell en Diebold, dit model wordt in Sectie 2.4.1 besproken. Cao and Wei (2000a) en Caballero et al. (2002) modelleren de temperatuur eveneens vanuit een econometrisch standpunt. Om deze gedetailleerde modellen voor de temperatuur te kunnen begrijpen, bespreken we eerst de fundamentele 7 8 Hoofdstuk 2. Modellering van de temperatuur econometrische modellen. Voorspellingen maken een belangrijk deel uit van de econometrische analyse, in deze sectie worden enkele voorspellingsmethoden besproken. Voor deze bespreking baseren we ons op Gujarati and Porter (2009) en Croux (2013). Om voorspellingen te doen, hebben we een model nodig, dit model moet de gegevens in de tijdreeks weergeven. Het proces dat we willen voorspellen moet dan voldoen aan het gevonden model. Aan de hand van regressie-analyse kunnen we een model opstellen dat de waarschijnlijke uitkomst van de tijdreeks in de nabije toekomst, gegeven de kennis van de meest recente uitkomsten, probeert te beschrijven. De econometrische modellen die gebruikt worden om de temperatuur te modelleren zijn voornamelijk autoregressieve (AR) modellen. Definitie 5. Een model is autoregressief wanneer het ´e´en of meerdere vorige waarden van de afhankelijke variabele (Y ) bevat naast de verklarende variabelen (X). Een voorbeeld van een autoregressief model is Yt = α + β · Xt + γ · Yt−1 + ut . Deze modellen worden ook vaak dynamische modellen genoemd, omdat ze een beeld vormen van het tijdspad van de afhankelijke variabele in relatie tot zijn voorbije waarden. We zullen nu de verschillende variaties iets meer in detail bespreken. Voor de meeste van deze variaties is stationariteit van het stochastisch proces Y vereist, we defini¨eren dit zoals in Croux (2013). Definitie 6. Zij Y = (Yt )t≥0 het stochastisch proces waaruit de tijdreeks verkregen wordt, dan is Yt (zwak) stationair als: • E[Yt ] gelijk is voor elke t; • Var(Yt ) gelijk is voor elke t; • Cov(Yt , Yt−k ) gelijk is voor elke t en voor elke k > 0. 2.1.1 MA-model In deze sectie defini¨eren we het voortschrijdend gemiddelde1 (MA) model. Veronderstel dat we een stationair stochastisch proces Y als volgt modelleren Yt = µ + ut + β1 · ut−1 , 1 In het Engels moving average. 2.1. Dynamische econometrische modellen 9 met µ een constante en u een ongecorreleerde stochastische storingsterm met gemiddelde 0 en constante variantie σ 2 (witte ruis). Hier is Y op tijdstip t gelijk aan een constante plus een voortschrijdend gemiddelde van de huidige en voorgaande storingstermen. In dit geval zeggen we dat Y een eerste orde voortschrijdend gemiddelde, MA(1), proces volgt. Wanneer Y echter voldoet aan Yt = µ + ut + β1 · ut−1 + β2 · ut−2 , dan is het een MA(2)-proces. Algemeen geldt dat een MA(q)-proces als volgt voorgesteld wordt Yt = µ + ut + β1 · ut−1 + β2 · ut−2 + · · · + βq · ut−q . Kortom, een MA-proces is gewoonweg een lineaire combinatie van witte ruis storingstermen. Tot slot geven we nog mee hoe men kan nagaan of een MA-proces een goed model is om een bepaalde tijdreeks te beschrijven. De autocorrelaties van een MA-proces gedragen zich op een bijzondere manier. Definitie 7. Voor een stationair stochastisch proces Y defini¨eren we de autocorrelatie van orde k als volgt Cov(Yt , Yt−k ) Cov(Yt , Yt−k ) ρk = Corr(Yt , Yt−k ) = p = . Var(Yt ) Var(Yt )Var(Yt−k ) De autocorrelaties kunnen geschat worden door PT (yt − y¯)(yt−k − y¯) . ρˆk = t=k+1 PT ¯)2 t=1 (yt − y De autocorrelaties van een MA(1)-proces worden in het bijzonder gegeven door ρ0 = 1, ρ1 = Corr(Yt , Yt−1 ) = β1 , 1 + β12 ρ2 = 0, ρ3 = 0, ... Dit kunnen we ook waarnemen in het correlogram dat hoort bij een MA(1)-proces. Definitie 8. Een grafiek die ρˆk plot tegenover k wordt een correlogram genoemd. 10 Hoofdstuk 2. Modellering van de temperatuur Op een correlogram kunnen we ook steeds twee lijnen waarnemen, deze corresponderen √ met de kritieke waarden van de teststatistiek ρˆk T om de nulhypothese ρk = 0 te testen voor een specifieke k. Op die manier weten we dat de tijdreeks afkomstig is van een MA(1)-proces wanneer enkel de waarden voor ρˆ0 en ρˆ1 ´e´en van deze lijnen overschrijden (zie Figuur 2.1). In het algemeen geldt voor de autocorrelaties van een MA(q)-proces Figuur 2.1: Correlogram van een MA(1)-proces. dat deze allen gelijk aan nul zijn voor lags groter dan q. Als het correlogram een sterke afname vertoont en niet significant wordt na lag q, kan men vermoeden dat de tijdreeks gegenereerd werd door een MA(q)-proces. We merken nog op dat we in alle volgende afbeeldingen van correlogrammen de autocorrelatie bij lag nul weglaten, want deze is toch steeds gelijk aan ´e´en. 2.1.2 AR-model Een stationair stochastisch proces Y = (Yt )t≥0 is autoregressief van orde ´e´en, AR(1), als het voldoet aan Yt = µ + α1 · Yt−1 + ut , met µ en α1 ongekende parameters en u de witte ruis stochastische storingsterm. De waarde van Y op tijdstip t is hier afhankelijk van zijn waarde op het vorige tijdstip enerzijds en van een stochastische term anderzijds. Dit model zegt met andere woorden dat de waarde van Y op tijdstip t, op een constante na, gewoon een proportie (= α1 ) van zijn waarde op tijdstip (t − 1) plus een willekeurige schok of storing op tijdstip t is. 2.1. Dynamische econometrische modellen 11 Wanneer we het volgende model beschouwen Yt = µ + α1 · Yt−1 + α2 · Yt−2 + ut , dan zeggen we dat Y een tweede orde autoregressief, AR(2), proces volgt. Dat is, de Y -waarde op tijdstip t hangt af van zijn waarde op de twee voorgaande tijdstippen, t − 1 en t − 2, en van een willekeurige storing op tijdstip t. In het algemeen hebben we: Y is een autoregressief proces van orde p, AR(p), wanneer Yt = µ + α1 · Yt−1 + α2 · Yt−2 + · · · + αp · Yt−p + ut geldt. Merk op dat in alle voorgaande modellen alleen de Y -waarden voorkomen, we hebben eventuele andere regressoren buiten beschouwing gelaten. Uit de autocorrelaties van een AR-proces kunnen we niet zo veel afleiden als voor een MA-proces. Zo worden de autocorrelaties voor een AR(1)-proces gegeven door ρ0 = 1, ρ1 = Corr(Yt , Yt−1 ) = α1 , ρ2 = α12 , ρ3 = α13 , ... Hier kunnen we niet veel uit afleiden. Voor AR-processen is het beter om naar de parti¨ele correlaties te kijken. Definitie 9. De parti¨ ele correlaties worden voor k = 0, 1, 2, . . . gegeven door πk = Corr(Yt , Yt−k |Yt−1 , . . . , Yt−k+1 ). Dit is gelijk aan de correlatie tussen Yt−k en de residuen van een regressie van Yt op de variabelen Yt−1 , . . . , Yt−k+1 . Specifiek voor een AR(1)-proces worden de parti¨ele autocorrelaties gegeven door π0 = 1, π1 = Corr(Yt , Yt−1 ) = α1 , π2 = 0, π3 = 0, ... Het partieel correlogram kan ons dus helpen om een AR(1)-proces te herkennen. 12 Hoofdstuk 2. Modellering van de temperatuur Definitie 10. Stel dat de parti¨ele autocorrelaties kunnen geschat worden door π ˆk . Een grafiek die π ˆk plot tegenover k wordt een partieel correlogram genoemd. Zoals bij het correlogram zijn er op het partieel correlogram ook twee lijnen aanwezig. √ Deze corresponderen met de kritieke waarden van de teststatistiek π ˆk T om de nulhypothese πk = 0 te testen voor een specifieke k. Op die manier weten we dat de tijdreeks afkomstig is van een AR(1)-proces wanneer enkel de waarden voor π ˆ0 en π ˆ1 ´e´en van deze lijnen overschrijden2 (zie Figuur 2.2). In het algemeen zijn de parti¨ele correlaties van lag groter dan p gelijk aan nul voor een AR(p)-proces. Figuur 2.2: Partieel correlogram van een AR(1)-proces. 2.1.3 ARMA-model Onderstel dat Y een stationair stochastisch proces is. Het is evident dat wanneer Y een ARMA-proces volgt, het kenmerken van zowel een AR- als een MA-proces heeft. Dus, Y volgt een ARMA(1,1)-proces, met ´e´en autoregressieve en ´e´en voortschrijdend gemiddelde term, als het kan geschreven worden als Yt = θ + α1 · Yt−1 + ut + β1 · ut−1 , waarbij θ een constante, α1 en β1 ongekende parameters en u de witte ruis stochastische storingsterm is. In het algemeen hebben we dat er in een ARMA(p,q)-proces enerzijds p 2 We merken op dat, gebruik makend van R, in het partieel correlogram de parti¨ele correlatie voor lag 0 niet weergegeven wordt. Uit het bovenstaande weten we dat steeds π0 = 1 geldt. 2.1. Dynamische econometrische modellen 13 autoregressieve en anderzijds q voortschrijdend gemiddelde termen zijn, Y voldoet dan aan Yt = θ + α1 · Yt−1 + · · · + αp · Yt−p + ut + β1 · ut−1 + · · · + βq · ut−q , met θ een constante. 2.1.4 ARIMA-model Zoals in Sectie 2.1.3 aangegeven werd, zijn ARMA-processen stationair. Soms geldt stationariteit voor een tijdreeks echter niet, maar is deze wel ge¨ıntegreerd van een bepaalde orde. Definitie 11. Y is ge¨ıntegreerd van orde ´ e´ en wanneer na het toepassen van de differentie-operator ∆ (∆Yt = Yt − Yt−1 ) een stationair proces verkregen wordt. Wanneer dus stationariteit verkregen wordt door de differentie-operator op de tijdreeks te laten inwerken, kunnen we een ARIMA-model specifi¨eren. Yt is een ARIMA(p,d,q)proces als ∆d Yt , met ∆d Yt = Yt −Yt−d , een ARMA(p,q)-proces is. De orde van integratie d is vaak nul of ´e´en. 2.1.5 SARIMA-model De algemene regel is dat wanneer er seizoenseffecten van orde s (vaak s = 12 of s = 4) aanwezig zijn, de seizoensgebonden differentie-operator ∆s gebruikt moet worden om stationariteit te verkrijgen. Vaak zal er echter nog enige seizoensgebondenheid in de correlatiestructuur van de tijdreeks achterblijven. Een sobere manier om zulke correlatiestructuren te modelleren is door gebruik te maken van SARMA-modellen. Een SARMA(0,1)(0,1)-model voor s = 12 wordt als volgt gespecifieerd Yt = c + (I − θ1 L)(I − Θ1 L12 )ut = c + (I − θ1 L)(ut − Θ1 ut−12 ) = c + ut − Θ1 ut−12 − θ1 ut−1 + θ1 Θ1 ut−13 , waarbij I de identiteitsoperator (I Yt = Yt ), L de lag-operator (L Yt = Yt−1 ) en u witte ruis is. Een SARMA(1,0)(0,1)-model voor s = 12 wordt als volgt gespecifieerd (I − φ1 L)Yt = c + (I − Θ1 L12 )ut ⇒Yt − φ1 Yt−1 = c + (ut − Θ1 ut−12 ) 14 Hoofdstuk 2. Modellering van de temperatuur ⇒Yt = c + φ1 Yt−1 + ut − Θ1 ut−12 , waarbij I de identiteitsoperator, L de lag-operator en u witte ruis is. In het algemeen geldt dat een stationair stochastisch proces Y een SARMA(p,q)(P ,Q)proces van orde s is als (I − φ1 L − . . . − φp Lp )(I − Φ1 Ls − . . . − Φp LsP )Yt = c + (I − θ1 L − . . . − θq Lq )(I − Θ1 Ls − . . . − ΘQ LsQ )ut , waarbij I de identiteitsoperator, L de lag-operator en u witte ruis is. Wanneer de tijdreeks niet stationair is, kunnen we –analoog aan de ARIMA-processen– de SARIMA-processen defini¨eren. De SARIMA-processen houden er rekening mee dat de tijdreeks een trend en seizoenseffecten kan vertonen, en op die manier dus niet stationair is. Immers, een lineaire trend kan verwijderd worden door de differentie-operator ∆ toe te passen en seizoenseffecten kunnen verwijderd worden door de seizoensgebonden differentie-operator ∆s te laten inwerken op de tijdreeks. Een SARIMA(p,d,q)(P ,D,Q)specificatie voor Yt komt dan ook overeen met een SARMA(p,q)(P ,Q)-model voor de variabele ∆d ∆D s Yt , het impliceert immers dat de differentie-operator d keer werd toegepast en de seizoensgebonden differentie-operator D keer werd gebruikt. 2.1.6 GARCH-model Een veralgemeend autoregressief conditioneel heteroscedastisch3 (GARCH) model is een model voor de met de tijd vari¨erende variantie van een bepaalde tijdreeks. Het eenvoudigste GARCH-model is het GARCH(1,1)-model, dat men als volgt kan schrijven 2 σt2 = α0 + α1 · u2t−1 + α2 · σt−1 . De conditionele variantie van u op tijdstip t hangt dus af van de gekwadrateerde storingsterm in de voorgaande tijdsperiode (u2t−1 ), maar ook van de conditionele variantie 2 in de voorgaande tijdsperiode (σt−1 ). Dit model kan natuurlijk veralgemeend worden tot een GARCH(p,q)-model dat p vorige waarden van de gekwadrateerde storingsterm en q vorige waarden van de conditionele variantie bevat, namelijk 2 2 . σt2 = α0 + α1 · u2t−1 + · · · + αp · u2t−p + β1 · σt−1 + · · · + βq · σt−q 3 In het Engels generalized autoregressive conditional heteroscedasticity. 2.2. CAR-model 2.2 15 CAR-model Het tijdscontinue autoregressieve (CAR) proces is een overgang tussen de econometrische aanpak en de tijdscontinue methode. Het proces vertrekt van het econometrische ARmodel en maakt het continu door differentialen te beschouwen. Een goede illustratie van deze link tussen econometrische en tijdscontinue modellen is het Ornstein-Uhlenbeckproces, dat hieronder besproken wordt, het kan immers beschouwd worden als het tijdscontinue analogon van het AR(1)-proces. In het algemeen zijn de CAR-processen de tijdscontinue analogi van de AR-processen. We volgen de definitie van Benth and Benth (2011), zij defini¨eren de temperatuur T (t) als volgt T (t) = Λ(t) + Υ(t) t ≥ 0, waarbij Λ(t) de functie is die het deterministische seizoensgemiddelde voorstelt en Υ een stochastisch proces is dat de fluctuaties rond het gemiddelde modelleert. Met andere woorden Υ(t) = T (t)−Λ(t) is de deseasonalized temperatuur. Uit empirisch onderzoek van de data blijkt dat Υ een autoregressieve structuur heeft op dagelijkse schaal. Dit motiveert het gebruik van een tijdscontinu autoregressief proces voor Υ. Voor p ≥ 1 wordt een CAR(p)-proces als volgt gedefinieerd. Definitie 12. Een tijdscontinu autoregressief proces van orde p wordt gedefinieerd als Υ(t) = e01 X(t), waarbij de notatie x0 de getransponeerde van x weergeeft. e01 is de canonische eerste eenheidsvector in Rp en X(t) voldoet aan (2.2), zie Sectie 2.3.1. 2.3 Tijdscontinue modellen Ge¨ınspireerd door Benth and Benth (2012), zullen we in deze sectie de tijdscontinue modellen bespreken. De processen die door deze modellen gedefinieerd worden, zijn immers het best geschikt wanneer we derivaten willen gaan prijzen. Daarnaast evolueert de temperatuur continu in de tijd, hierdoor is het zeer logisch om een tijdscontinu stochastisch proces te gebruiken om de temperatuurschommelingen te modelleren ook al zullen de gegevens vaak op dagelijkse basis verkregen worden (bv. gemiddelde dagtemperatuur). In tijdscontinue modellering wordt het proces verkregen als een oplossing van een stochastische differentiaalvergelijking. Omdat wij ge¨ınteresseerd zijn in het verkrijgen van 16 Hoofdstuk 2. Modellering van de temperatuur een temperatuurproces zullen we modellen beschouwen die de mean-reverting eigenschap bezitten, cf. Alaton et al. (2002). Deze specificatie is ingegeven door het volgende kenmerk van de temperatuur: het is niet mogelijk dat de temperatuur dag na dag blijft stijgen (of dalen) gedurende een lange periode. Dit wil zeggen dat de temperatuur altijd zal schommelen rond de gemiddelde waarde en na een stijging of daling ook altijd zal terugkeren naar dit gemiddelde. Hierdoor zal ons model enkel afwijkingen van de gemiddelde waarde over korte tijdsintervallen mogen toelaten, en dit wordt uitgedrukt door de mean-reverting eigenschap. We bespreken de twee meest gebruikte modellen om het temperatuurproces continu te beschrijven, het Ornstein-Uhlenbeck-model en het Hull-White-model. De literatuur is niet altijd even duidelijk over het onderscheid tussen beide modellen, wij gebruiken de volgende definities. 2.3.1 Ornstein-Uhlenbeck Aan de basis van het Ornstein-Uhlenbeck-model ligt, volgens Smith (2010), de volgende stochastische differentiaalvergelijking dX(t) = λ · (µ − X(t))dt + σdW (t), (2.1) met positieve parameters λ, µ en σ en {W (t), t ≥ 0} een Brownse beweging. Het komt vaak voor dat λ = 0 wordt gekozen, cf. Valdivieso et al. (2009). Het stochastisch proces X = {X(t), t ≥ 0} wordt verkregen door de stochastisch differentiaalvergelijking op te lossen met behulp van een integrerende factor. We vermenigvuldigen beide leden van dX(t) = λ · (µ − X(t))dt + σdW (t) ⇔ dX(t) + λX(t)dt = µλdt + σdW (t), met eλt om te komen tot eλt dX(t) + λeλt X(t)dt = eλt · (µλdt + σdW (t)) d(eλt X(t)) = µλeλt dt + σeλt dW (t). ⇔ De laatste vergelijking is nu onmiddellijk integreerbaar, we integreren over [0, T ] Z T Z T Z T λt λt d(e X(t)) = µλe dt + σeλt dW (t) 0 ⇔ ⇔ λT e X(T ) − e 0 λ·0 0 λt e X(0) = µλ λ T T Z σeλt dW (t) + 0 eλT X(T ) − X(0) = µ · (eλT − eλ·0 ) + 0 Z 0 T σeλt dW (t) 2.3. Tijdscontinue modellen 17 e X(T ) = X(0) + µ · (e λT Z T σeλt dW (t) 0 Z T −λT −λT −λT eλt dW (t), ⇔ X(T ) = e X(0) + µ · (1 − e ) + σe ⇔ λT − 1) + 0 als X(T ) van deze vorm is, volgt X een Ornstein-Uhlenbeck-proces. In Benth and Benth (2011) wordt gebruik gemaakt van het p-dimensionaal OrnsteinUhlenbeck-proces om een CAR(p)-proces te defini¨eren. Definieer de p × p-matrix A als volgt 0 1 0 ··· 0 0 0 1 · · · 0 . .. .. .. .. A = .. . . . . , 0 0 0 ··· 1 −ap −ap−1 −ap−2 · · · −a1 waarbij ai , i = 1, . . . , p, positieve constanten zijn. Voor een ´e´en-dimensionale Brownse beweging {W (t), t ≥ 0} wordt het p-dimensionaal Ornstein-Uhlenbeck-proces X gedefinieerd als dX(t) = AX(t)dt + σep dW (t), (2.2) 0 0 . waarbij σ > 0 constant en ep = .. , de canonische p-de eenheidsvector in Rp . De 0 1 oplossing {X(t), t ≥ 0} van (2.2) voldoet aan Z t At X(t) = e X(0) + eA(t−u) σep dW (u). 0 Dit kan ingezien worden door dezelfde redenering te volgen als voor (2.1) en −λ = A, µ = 0 en σ = ep σ te kiezen. 2.3.2 Hull-White In het Hull-White-model wordt het stochastisch proces verkregen door het oplossen van volgende stochastische differentiaalvergelijking (cf. Shreve (2004)) dX(t) = (α(t) − β(t)X(t))dt + σ(t)dW (t), (2.3) 18 Hoofdstuk 2. Modellering van de temperatuur met α(t), β(t) en σ(t) niet stochastische, positieve functies van de tijdsveranderlijke t en {W (t), t ≥ 0} een Brownse beweging. Om de vergelijking dX(t) = (α(t) − β(t)X(t))dt + σ(t)dW (t) ⇔ dX(t) + β(t)X(t)dt = α(t)dt + σ(t)dW (t), op te lossen, maken we opnieuw gebruik van een integrerende factor en vermenigvuldigen hiertoe beide leden met e Rt e 0 β(s)ds 0 β(s)ds dX(t) + β(t)e Rt ⇔ Rt d(e 0 β(s)ds Rt 0 β(s)ds X(t)dt = e Rt X(t)) = α(t)e 0 β(s)ds Rt 0 β(s)ds · (α(t)dt + σ(t)dW (t)) Rt dt + σ(t)e 0 β(s)ds dW (t). Dit resultaat is onmiddellijk integreerbaar, integratie over [0, T ] levert ⇔ ⇔ ⇔ 0 β(s)ds R0 Z X(T ) − e 0 β(s)ds T Rt Z β(s)ds T Rt X(0) = α(t)e 0 dt + σ(t)e 0 β(s)ds dW (t) 0 0 Z T Z T Rt Rt RT σ(t)e 0 β(s)ds dW (t) α(t)e 0 β(s)ds dt + e 0 β(s)ds X(T ) − X(0) = 0 0 Z T Z T Rt Rt RT σ(t)e 0 β(s)ds dW (t) α(t)e 0 β(s)ds dt + e 0 β(s)ds X(T ) = X(0) + 0 0 Z T Z T Rt Rt RT β(s)ds β(s)ds − 0 β(s)ds σ(t)e 0 dW (t) α(t)e 0 dt + X(T ) = e X(0) + 0 0 Z T Z T R R RT − 0T β(s)ds − tT β(s)ds X(T ) = e X(0) + α(t)e dt + σ(t)e− t β(s)ds dW (t), e ⇔ RT 0 0 X(T ) voldoet aan (2.3), het is een Hull-White-proces. 2.4 Praktijk In deze sectie zullen we in het kort enkele modellen bespreken die in de praktijk gebruikt worden om de temperatuur te modelleren. Het Alaton-model dat we hieronder uitgebreider zullen bespreken, is een model dat de temperatuur als een tijdscontinu proces beschouwt. Het Benth-model gaat ditzelfde proces discretiseren, maar is nog steeds te plaatsen onder de continue modellen. De keuze van deze twee modellen werd ingegeven door Schiller et al. (2012). Om te beginnen bespreken we een model dat te plaatsen is bij de modellen die een discreet proces gebruiken om de dagtemperatuur te modelleren, namelijk het model voorgesteld door Campbell en Diebold. 2.4. Praktijk 2.4.1 19 Het model van Campbell en Diebold Het is niet noodzakelijk om over een structureel model te beschikken wanneer men succesvol wil voorspellen. In Campbell and Diebold (2005) vertrekt men van deze gedachte en vraagt men zich af of het mogelijk is om de temperatuur te modelleren en te voorspellen via een niet-structurele tijdreeksbenadering. Campbell en Diebold modelleren en voorspellen direct de gemiddelde dagtemperatuur4 (Temp) gemeten in graden Fahrenheit en doorlopen hiervoor de stappen om een model te bouwen. Eerst en vooral analyseren ze de gegevens, het is immers nuttig om een algemene indruk over de temperatuurdata te krijgen. De grafiek van de tijdreeks vertoont een sterke aanwezigheid van seizoenseffecten in de temperatuur, de temperatuur beweegt herhaaldelijk en regelmatig door perioden van hoge temperatuur (zomer) en lage temperatuur (winter). Tot nu toe blijkt dat een seizoensgebonden component belangrijk zal zijn in ieder tijdreeksmodel dat de temperatuur beschrijft, omdat de temperatuur een uitgesproken seizoensgebonden variatie bezit. Campbell en Diebold verkiezen het gebruik van een Fourierreeks van lage orde om deze seizoensgevoeligheid te modelleren, omdat dit twee voordelen heeft. Enerzijds produceert het een glad seizoenspatroon, dit is in overeenstemming met de fundamentele intu¨ıtie dat het doorlopen van de verschillende seizoenen eerder gradueel in plaats van discontinu gebeurt. Anderzijds bevordert het parsimony (de drang om verschijnselen te verklaren met behulp van minder parameters), wat de complexiteit van het model verlaagt en dus de numerieke stabiliteit van schattingen verhoogt. Men kan vermoeden dat naast deze seizoenseffecten ook niet-seizoensgebonden factoren kunnen meespelen in de dynamiek van de temperatuur. Een dergelijke factor is de trend, deze kan relevant zijn maar is hier waarschijnlijk gering, gezien de korte (40 jaar) observatieperiode die Campbell en Diebold hanteren. Zij kiezen dus voor een eenvoudige polynomiale deterministische trend van lage orde. Nog zo’n factor is de cyclus waarmee ze elk ander persistent verschijnsel, afgezien van de seizoenseffecten en de trend, bedoelen. Deze cyclische dynamiek leggen ze vast door gebruik te maken van autoregressieve lags. Tot nu toe lag de focus op de dynamiek van het conditioneel gemiddelde, met bijdragen van de trend, de seizoenseffecten en de cyclische component. Campbell en Diebold nemen echter ook de dynamiek van de conditionele variantie (volatiliteit) in aanmerking, met 4 Zoals in Benth and Benth (2012), nemen we als referentie voor de gemiddelde dagtemperatuur het gemiddelde van de maximum- en minimumtemperatuur, zoals in Definitie 1. Wij zullen in deze sectie dus “temperatuur”gebruiken waar Campbell en Diebold “gemiddelde dagtemperatuur”gebruiken. 20 Hoofdstuk 2. Modellering van de temperatuur bijdragen uit zowel de seizoens- als de cyclische component. Ze benaderen de seizoensgebonden volatiliteit aan de hand van een Fourierreeks en de cyclische volatiliteit aan de hand van een veralgemeend autoregressief conditioneel heteroscedastisch (GARCH) proces, zie Sectie 2.1.6. Na het analyseren van de gegevens stellen ze een model voor, dat de besproken kenmerken bevat. Ze schatten het volgende model voor de temperatuur op dag t Tempt = Trendt + Seasonalt + L X ρt−l Tempt−l + σt εt , (2.4) l=1 met Trendt = M X βm tm , m=0 P X d(t) d(t) Seasonalt = σc,p cos 2πp + σs,p sin 2πp , 365 365 p=1 X Q S K X X d(t) d(t) 2 2 2 βs σt−s , αk (σt−k εt−k ) + γc,q cos 2πq + γs,q sin 2πq + σt = 365 365 s=1 q=1 k=1 εt ∼ iid(0, 1), waarbij d(t) een repeterende trapfunctie is die door 1, . . . , 365 loopt (laat in alle schrikkeljaren 29 februari vallen). Campbell en Diebold kiezen L = 25, M = 1, P = 3, Q = 3, K = 1 en S = 1, op basis van zowel het Akaike-informatiecriterium (AIC) als het Schwarz-informatiecriterium (SIC). Men wil de AIC- en SIC-waarde zo klein mogelijk houden, met p AIC = ln(MSE) + 2 , T p ln(T ) , SIC = ln(MSE) + T waarbij p het aantal parameters in het model aangeeft en MSE staat voor de Mean T 1X 2 (t−1) Squared Error uˆ met uˆt de one-step ahead forecast errors uˆt = yt − yˆt T t=1 t (cf. Croux (2013)). Op die manier penaliseren het AIC en het SIC (zelfs meer dan het AIC) voor de complexiteit van het model. Het regressiemodel wordt geschat aan de hand van de quasi-maximum-likelihoodmethode, volgens Bollerslev and Wooldridge (1992) gaat dit als volgt. Zij {(yt , zt ) : t = 1, 2, . . .} een rij van observeerbare stochastische vectoren met yt en zt respectievelijk een K × 1 en een L × 1 matrix. Stel dat de vector yt de endogene en de vector zt de exogene variabelen 2.4. Praktijk 21 bevat. Zij xt = (zt , yt−1 , zt−1 , . . . , y1 , z1 ) de vooraf bepaalde variabelen. Het conditionele gemiddelde en de variantie zijn gezamenlijk geparametriseerd door een eindigdimensionale vector θ {µt (xt , θ) : θ ∈ Θ}, {Ωt (xt , θ) : θ ∈ Θ}, waarbij Θ een deelverzameling van Rp is en waarbij µt en Ωt gekende functies van xt en θ zijn. Voor observatie t geldt er dat de quasi-conditionele log-likelihood gelijk is aan (op een constante na) 1 1 lt (θ; yt , xt ) = − ln |Ωt (xt , θ)| − (yt − µt (xt , θ))0 Ω−1 t (xt , θ)(yt − µt (xt , θ)). 2 2 Stel nu εt (yt , xt , θ) ≡ yt − µt (xt , θ) de K × 1 residuele vector. Door het achterwege laten van de xt - en yt -afhankelijkheid in εt en Ωt , bekomen we de compactere notatie 1 1 lt (θ) = − ln |Ωt (θ)| − εt (θ)0 Ω−1 t (θ)εt (θ). 2 2 De quasi-maximum-likelihoodschatter (QMLE) θˆT wordt verkregen door het maximaliseren van de quasi-log-likelihoodfunctie LT (θ) = T X lt (θ). t=1 Na de schatting van de parameters in het model, moet het model nog gevalideerd worden. Campbell en Diebold concluderen dat voor het model dat de dynamiek van het conditionele gemiddelde beschrijft (vergelijking (2.4) samen met de vergelijkingen voor Trendt en Seasonalt ) alle termen statistisch significant zijn. Zowel de trend, de seizoensgevoeligheid als de cyclische persistentie zijn significant. De residuen voldoen eveneens aan de nodige voorwaarden. Daarbovenop vertoont dit deel van het model een determinatieco¨effici¨ent5 R2 die groter is dan 90%, wat betekent dat meer dan 90% van de variatie van de temperatuur verklaard wordt door het model. Het model is dus goed aangepast aan de data. Uit bepaalde grafieken is de seizoensgebonden conditionele heteroscedasticiteit (de conditionele variantie is niet constant) in de temperatuur af te leiden. Het model voor de volatiliteit (de formule voor σt2 ) werd ontworpen om deze conditionele heteroscedasticiteit te benaderen. Uit de schatting van dit model blijkt dat beide termen, zowel 5 De determinatieco¨effici¨ent R2 is een maat die informatie geeft over de mate waarin een model de werkelijke data benadert. 22 Hoofdstuk 2. Modellering van de temperatuur de seizoensgebonden (Fourier) als de cyclische (GARCH) component, significant zijn. Campbell en Diebold vinden dat de component van de seizoensgebonden volatiliteit (Fourier) relatief gezien belangrijker is, hij is significant en tamelijk groot. De nietseizoensgebonden GARCH-volatiliteit heeft een kleiner effect, de persistentie wordt gemeten door de GARCH-parameter β. Om af te sluiten hebben Campbell en Diebold het correlogram van de gekwadrateerde gestandaardiseerde residuen6 geplot, waaruit ze konden concluderen dat er geen significant verschil was met het gedrag van witte ruis. Met andere woorden, het door hen gekozen model is voldoende om de data te beschrijven. 2.4.2 Het Alaton-model Volgens Schiller et al. (2012), die zich tevens baseren op Alaton et al. (2002), past het Alaton-model de parameters van het Hull-White-model (zie Sectie 2.3.2) aan zodanig dat ze de kenmerken van de temperatuur weerspiegelen. Zij S(t) een deterministische functie die de seizoensgevoeligheid modelleert, dan kan deze als volgt gedefinieerd worden S(t) = A + Bt + C sin(ωt + φ). (2.5) Deze keuze is afkomstig van het feit dat Alaton et al. (2002) in hun temperatuurgegevens een sterke seizoensgebonden variatie terugvinden, die te modelleren is aan de hand van een sinusfunctie sin(ωt + φ). Hier staat t voor de tijd uitgedrukt in dagen en voor de pe2π riode geldt ω = , omdat het seizoensgebonden karakter van de temperatuur jaarlijks 365 terugkeert en de temperatuur op dagelijkse basis gemodelleerd wordt (we verwaarlozen de schrikkeldagen, in de praktijk worden ze uit de gegeven data verwijderd). De fase φ ∈ [−π, π] geeft weer dat de jaarlijkse minimum- en maximumtemperatuur niet noodzakelijk, respectievelijk, op 1 januari en 1 juli voorkomen7 . Daarnaast vertoont hun data ook een positieve trend, d.w.z. dat de temperatuur jaarlijks stijgt. Deze stijging is heel klein, waardoor ze veronderstellen dat ze de trend lineair kunnen voorstellen (A + Bt). De parameter C definieert de amplitude van het verschil tussen de jaarlijkse minimale en maximale gemiddelde dagtemperatuur. De parameters A, B, C en φ moeten zodanig gekozen worden dat de curve de data op een correcte manier weergeeft. Omdat de temperatuur niet deterministisch is, moeten we ook rekening houden met eventuele ruis. Alaton et al. (2002) kiezen ervoor om hiertoe de Brownse beweging 6 De gestandaardiseerde residuen worden weergegeven door Tt −Tbt σ ˆt , waarbij Tbt de gefitte waarde van de temperatuur is. 7 Het interval [−π, π] is het enige relevante interval dat men moet doorlopen, omdat een cos-functie eigenlijk een sin-functie met een fase is. 2.4. Praktijk 23 {W (t), t ≥ 0} te gebruiken. Op die manier wordt de temperatuur op tijdstip t, T (t), gemodelleerd aan de hand van de volgende specifieke stochastische differentiaalvergelijking dS(t) dT (t) = a · (S(t) − T (t)) + dt + σ(t)dW (t), dt t ≥ 0, (2.6) met a de snelheid van mean-reversion, σ(t) de seizoensgevoeligheid van de dagelijkse temperatuurverandering van de residuen en {W (t), t ≥ 0} een Brownse beweging. De oplossing van deze stochastische differentiaalvergelijking kan bekomen worden zoals in Sectie 2.3.2 met X(t) = T (t) − S(t), α(t) = 0 en β(t) = a, dit levert Z t −at T (t) = S(t) + e · (T (0) − S(0)) + e−a(t−s) σ(s)dW (s). (2.7) 0 Om een bruikbaar model te verkrijgen, moeten de parameters geschat worden, zodanig dat ze de data goed weergeven. Om de numerieke waarden van de constanten in (2.5) te vinden, passen Alaton et al. (2002) de volgende functie aan de temperatuurdata aan Yt = a1 + a2 t + a3 sin(ωt) + a4 cos(ωt), (2.8) en dit volgens de methode van de kleinste kwadraten. Dit betekent dat we de parametervector ξ = (a1 a2 a3 a4 ) moeten vinden die een oplossing is van min ||Y − X||2 , ξ waarbij Y de vector met elementen (2.8) is en X de vector die de data bevat. De constanten in (2.5) worden dan als volgt verkregen S(t) = Yt ⇔ A + Bt + C sin(ωt + φ) = a1 + a2 t + a3 sin(ωt) + a4 cos(ωt) ⇒ A = a1 , B = a2 , C sin(ωt + φ) = a3 sin(ωt) + a4 cos(ωt). Uit de laatste gelijkheid kunnen we de schatting voor C en φ als volgt bepalen C sin(ωt + φ) = a3 sin(ωt) + a4 cos(ωt) ⇔ C[sin(ωt) cos(φ) + cos(ωt) sin(φ)] = a3 sin(ωt) + a4 cos(ωt) 24 Hoofdstuk 2. Modellering van de temperatuur ( ⇒ C sin(ωt) cos(φ) = a3 sin(ωt) C cos(ωt) sin(φ) = a4 cos(ωt) ( C cos(φ) = a3 ⇒ . C sin(φ) = a4 We merken op dat het koppel (φ, C) eigenlijk het koppel poolco¨ordinaten van P (a3 , a4 ) is, er geldt dus a tan(φ) = 4 q a3 , C = a2 + a2 4 3 met C ≥ 0 en φ ∈ [−π, π]. Hieruit volgt dat de constanten in (2.5) gegeven worden door A = a1 , B = a2 , q C = a23 + a24 , a4 Bgtan als a3 > 0 a3 a4 Bgtan + π als a3 < 0, a4 ≥ 0 a3 a4 . φ = − π als a3 < 0, a4 < 0 Bgtan a3 π als a3 = 0, a4 > 0 2π − als a3 = 0, a4 < 0 2 Specifiek voor het Alaton-model wordt er verondersteld dat de variantie σ 2 (t) maande2 , gebruiken lijks ongeveer constant blijft. Als schatter van de maandelijkse variantie, σm we de kwadratische variatie van T (t), d.w.z. 2 ˆ¯m,y σ Nm −1 1 X = (T (t + 1)m,y − T (t)m,y )2 , Nm t=0 2 σ ˆm Ny 1 X 2 ˆ¯ , = σ Ny y=1 m,y met Nm het aantal dagen in maand m, T (t)m,y de temperatuur op dag t in maand m in jaar y en Ny het aantal gebruikte jaren. Nu moet enkel nog de mean-reversion parameter a geschat worden, we volgen de aanpak van Alaton et al. (2002) die de martingale estimation functions methode gebruiken. 2.4. Praktijk 25 Normaal wordt parameterinferentie gebaseerd op de likelihoodfunctie L, de maximumlikelihoodschatter heeft immers meestal de nodige eigenschappen. Inferentie voor discrete ˜ van de continue likelihoodfuncobservaties kan gebaseerd worden op een benadering, L, ˜ verkregen wordt, is niet consistent wanneer de tijd tie. Echter, de schatter die via L tussen de observaties een ondergrens heeft die strikt groter is dan nul. Daarnaast kan de schatter vertekend zijn, wanneer de tijd tussen de observaties te groot is. Deze gebreken van de schatter worden vermeden door de constructie van een martingale estimation ˜ zie ook Bibby and Sørensen (1995). Het is toepasselijk om deze mefunction van L, thode hier te gebruiken, want de tijd tussen de observaties van de temperatuur (meestal een dag) is duidelijk strikt groter dan nul. De methode werkt als volgt. Een effici¨ente schatter a ˆn voor a (die gebaseerd is op observaties die over n dagen verzameld werden), voldoet aan de volgende vergelijking Gn (ˆ an ) = 0, (2.9) waarin Gn als volgt gedefinieerd wordt Gn (a) = n X b0 (T (i − 1); a) i=1 2 σi−1 (T (i) − E[T (i)|T (i − 1)]) , met b0 (T (t); a) de afgeleide naar a van de driftterm b(T (t); a) = a · (S(t) − T (t)) + Gn is een schattingsfunctie die een martingaal is, immers (2.10) dS(t) . dt E[Gn+1 (a) | n, n − 1, n − 2, . . .] # " n+1 X b0 (T (i − 1); a) =E (T (i) − E[T (i)|T (i − 1)]) n 2 σ i−1 " i=1 n X b0 (T (i − 1); a) =E (T (i) − E[T (i)|T (i − 1)]) 2 σ i−1 i=1 b0 (T (n); a) + (T (n + 1) − E[T (n + 1)|T (n)]) n 2 σn b0 (T (n); a) n (T (n + 1) − E[T (n + 1)|T (n)]) = E Gn (a) + σn2 b0 (T (n); a) = Gn (a) + (E[T (n + 1) | n] − E[T (n + 1)|T (n)]) σn2 = Gn (a). Om (2.9) op te lossen, is het voldoende dat we alle termen E[T (i)|T (i − 1)] in de definitie van Gn bepalen. Via (2.7), hebben we voor t ≥ s −a(t−s) T (t) = S(t) + e Z · (T (s) − S(s)) + s t e−a(t−τ ) σ(τ )dW (τ ), 26 Hoofdstuk 2. Modellering van de temperatuur dit levert voor t = i en s = i − 1 Z i e−a(i−τ ) σ(τ )dW (τ ) · (T (i − 1) − S(i − 1)) + i−1 Z i e−a(i−τ ) σ(τ )dW (τ ). = S(i) + e−a · (T (i − 1) − S(i − 1)) + T (i) = S(i) + e −a(i−(i−1)) i−1 Via de eigenschappen van de conditionele verwachtingswaarde, krijgen we het volgende E[T (i)|T (i − 1)] Z −a = E S(i) + e · (T (i − 1) − S(i − 1)) + i e −a(i−τ ) i−1 σ(τ )dW (τ )T (i − 1) = E[S(i)|T (i − 1)] + E[e−a · (T (i − 1) − S(i − 1))|T (i − 1)] Z i −a(i−τ ) +E e σ(τ )dW (τ )T (i − 1) i−1 Z i−1 −a = E[S(i)|T (i − 1)] + e · (T (i − 1) − S(i − 1)) + e−a(i−1−τ ) σ(τ )dW (τ ) i−1 −a = S(i) + e (T (i − 1) − S(i − 1)). De voorlaatste gelijkheid geldt door de martingaaleigenschap van de stochastische integraal en de laatste door het deterministisch karakter van S(t). We kunnen dit nu invullen in (2.10), er komt Gn (a) = n X S(i − 1) − T (i − 1) i=1 2 σi−1 (T (i) − (T (i − 1) − S(i − 1))e−a − S(i)), nu moeten we hier enkel nog de nulwaarde van zoeken om een schatter voor a te vinden. S(i − 1) − T (i − 1) We zullen voor de eenvoud van notatie Yi−1 = noteren. Stel Gn (a) = 2 σi−1 0, dit levert n X Yi−1 (T (i) − (T (i − 1) − S(i − 1))e−a − S(i)) = 0 i=1 ⇔ n X −a Yi−1 (T (i − 1) − S(i − 1))e i=1 = n X Yi−1 (T (i) − S(i)) i=1 Pn Yi−1 (T (i) − S(i)) ⇔ e−a = Pn i=1 Yi−1 (T (i − 1) − S(i − 1)) i=1 Pn Yi−1 (T (i) − S(i)) i=1 ⇔ −a = ln Pn . i=1 Yi−1 (T (i − 1) − S(i − 1)) We verkrijgen dus de volgende schatter Pn Yi−1 (T (i) − S(i)) i=1 a ˆn = − ln Pn i=1 Yi−1 (T (i − 1) − S(i − 1)) 2.4. Praktijk 27 Pn S(i−1)−T (i−1) (T (i) 2 i=1 σi−1 = − ln Pn S(i−1)−T (i−1) (T (i 2 i=1 σi−1 − S(i)) − 1) − S(i − 1)) . Na de schatting van de parameters, is het mogelijk om het proces te gaan simuleren (bv. via Monte Carlo-simulatie) en zo de distributie van de temperatuur te voorspellen. 2.4.3 Het Benth-model Volgens Schiller et al. (2012), die zich tevens baseren op Benth and Benth (2007), maakt het Benth-model ook gebruik van de stochastische differentiaalvergelijking (2.6). Het verschil met het Alaton-model zit hem enerzijds in de specificatie van de parameters S(t) en σ(t). De deterministische functie die de seizoensgevoeligheid en de trend specifieert, S(t), wordt hier gedefinieerd als een afgeknotte Fourierreeks met lineaire trend X I1 J1 X 2iπ(t − fi ) 2jπ(t − gj ) S(t) = b + ct + a0 + ai sin + . bj cos 365 365 i=1 j=1 De variantie, σ 2 (t), wordt in het Benth-model als volgt gespecifieerd X I2 J2 X 2iπt 2jπt 2 σ (t) = d + + . ci sin dj cos 365 365 i=1 j=1 Merk op dat ook in dit model de schrikkeldagen verwaarloosd worden (in de praktijk worden de schrikkeldagen uit de gegevens verwijderd). Daarnaast merken we eveneens op dat we σ 2 (t) als een afgeknotte Fourierreeks defini¨eren en niet σ(t) zelf, dit is het meest geschikt voor de data-analyse. Daarnaast is het ook beter vanuit het oogpunt om derivaten te prijzen, omdat σ 2 (t) vaker voorkomt dan σ(t). In Benth and Benth (2007) worden I1 = 1 en J1 = 1 gesteld om de functie S(t) aan de data aan te passen. De numerieke waarden die hierdoor voor a0 , a1 en b1 verkregen worden, zijn voldoende om de seizoensgebonden variatie van de temperatuur te beschrijven. Ze vinden eveneens dat het voldoende is om I2 = 4 en J2 = 4 te kiezen wanneer de variantie van de temperatuur bepaald moet worden, dit levert immers gepaste numerieke waarden voor d, ci en di , i = 1, . . . , 4. Zoals bij het Alaton-model, krijgen we de volgende oplossing van de stochastische differentiaalvergelijking T (t) = S(t) + e −at Z · (T (0) − S(0)) + t e−a(t−s) σ(s)dW (s). 0 Hieruit kunnen we afleiden dat de temperatuur op elk tijdstip t normaal verdeeld is, maar terugkeert naar het gemiddelde S(t). De snelheid van mean-reversion is a, en de 28 Hoofdstuk 2. Modellering van de temperatuur variantie van de temperatuur varieert seizoensgebonden met σ(t). Naast de specificaties van de parameters is er nog een ander verschil met het Alatonmodel, in het Benth-model wordt het temperatuurproces beschouwd als een discreet proces. Aangezien er wordt vertrokken van dezelfde stochastische differentiaalvergelijking (2.6), kan het discrete proces verkregen worden door discretisatie van (2.7). Dit gaat als volgt, beschouw T (t + 1) en T (t) −a(t+1) T (t + 1) = S(t + 1) + e Z t+1 e−a(t+1−s) σ(s)dW (s) · (T (0) − S(0)) + 0 Z t+1 e−a(t−s) σ(s)dW (s), · (T (0) − S(0)) + e 0 Z t T (t) = S(t) + e−at · (T (0) − S(0)) + e−a(t−s) σ(s)dW (s), −a −at −a = S(t + 1) + e e 0 dit impliceert T (t + 1) − T (t) = S(t + 1) − S(t) + (e−a − 1)e−at · (T (0) − S(0)) Z t Z t+1 −a −a(t−s) −a +(e − 1) e σ(s)dW (s) + e e−a(t−s) σ(s)dW (s). 0 t Uit (2.7) volgt dat (e −a − 1)e −at −a · (T (0) − S(0)) + (e Z − 1) t e−a(t−s) σ(s)dW (s) = (e−a − 1) · (T (t) − S(t)), 0 hierdoor kunnen we T (t + 1) − T (t) op de volgende manier herschrijven Z t+1 −a −a T (t + 1) − T (t) = S(t + 1) − S(t) + (e − 1) · (T (t) − S(t)) + e e−a(t−s) σ(s)dW (s). t Noteren we nu ∆X(t) = X(t + 1) − X(t), dan komt er −a ∆T (t) = ∆S(t) − (1 − e ) · (T (t) − S(t)) + e −a Z t+1 e−a(t−s) σ(s)dW (s). t Zoals in Einmahl (2012) kan men stochastische integralen via een soort Riemann-som benaderen Z T f (s)dW (s) ≈ 0 n X f (ti−1,n ) · (W (ti,n ) − W (ti−1,n )), i=1 Z t+1 T ·i waarbij ti,n = . Wij willen de integraal e−a(t−s) σ(s)dW (s) benaderen, uit het n t bovenstaande volgt dat dit op de volgende manier kan Z t+1 e−a(t−s) σ(s)dW (s) ≈ e−a(t−t) σ(t) · (W (t + 1) − W (t)) = σ(t)∆W (t). t 2.4. Praktijk 29 Via de benadering voor de stochastische integraal vinden we dus ∆T (t) ≈ ∆S(t) − (1 − e−a ) · (T (t) − S(t)) + e−a σ(t)∆W (t). Wanneer we deze uitdrukking vereenvoudigen, komt er T (t + 1) − T (t) ≈ S(t + 1) − S(t) − (1 − e−a ) · (T (t) − S(t)) + e−a σ(t)∆W (t) ⇔ T (t + 1) − S(t + 1) ≈ T (t) − S(t) − (1 − e−a ) · (T (t) − S(t)) + e−a σ(t)∆W (t) ≈ e−a · (T (t) − S(t)) + e−a σ(t)∆W (t). Stellen we nu Te(t) = T (t) − S(t), α = e−a en σ ˜ (t) = ασ(t), dan krijgen we de onderstaande tijdreeks Te(t + 1) = αTe(t) + σ ˜ (t)εt , waarbij εt i.i.d. standaard normaal verdeeld is. Dat εt normaal verdeeld moet zijn volgt uit ∆W (t) = W (t + 1) − W (t) en de definitie van een Brownse beweging, deze vereist dat voor 0 ≤ s < t < ∞ geldt dat W (t) − W (s) ∼ N (0, t − s). Om deze tijdreeks aan de data aan te passen, gaat men te werk in verschillende stappen. Empirische studies tonen aan dat data van de dagtemperatuur een significante lineaire trend vertonen, dus moet er eerst gecontroleerd worden of er in de gegeven data ook zo een trend aanwezig is. Vervolgens kunnen we S(t) numeriek bepalen via de methode van de kleinste kwadraten. Hierboven hebben we reeds vermeld dat de keuze van I1 = 1 en J1 = 1 voldoende is om een goede weergave van de data te krijgen. Daarna, voeren we een lineaire regressie uit die de temperatuur van vandaag tegenover de temperatuur van de voorgaande dagen plaatst. Ten slotte kunnen we de dagelijkse variantie van de regressie residuen schatten om de functie σ ˜ 2 (t) te bekomen. Merk op dat we hierboven al vermeld hebben dat de keuze van I2 = 4 en J2 = 4 voor de data in Benth and Benth (2007) voldoende was om σ 2 (t) te bepalen. De enige parameter die dus nog moet geschat worden in dit discrete proces is a, en dit kan eveneens met behulp van lineaire regressie. 2.4.4 Econometrische modellen voor Belgische temperatuurdata In deze sectie zullen we zelf twee modellen opstellen voor Belgische temperatuurdata. Het gaat hier over de zogenaamde dynamische econometrische modellen, die kunnen gespecifieerd worden via de analyse van een tijdreeks. We zullen enkele modellen specifi¨eren, schatten en valideren. We volgen de aanpak die aangeleerd werd in het vak Statistiek voor Actuari¨ele Wiskunde (cf. Croux (2013)). 30 Hoofdstuk 2. Modellering van de temperatuur De temperatuurgegevens waarmee we werken, werden geleverd door het KMI8 . Het zijn gegevens over de maximale en minimale temperatuur per dag gemeten in Ukkel van 01/01/1981 tot en met 31/12/2010. Om een model te specifi¨eren voor de temperatuur (zie Definitie 1), moeten we eerst de gegevens naar de gewenste vorm transformeren. We hebben informatie over de maximale (Timax ) en minimale (Timin ) temperatuur per dag gekregen en kunnen dus de temperatuur (Ti ) op dag i als volgt berekenen Ti = Timax + Timin . 2 Om op een goede manier een model te kunnen specifi¨eren, moet de tijdreeks regelmatig gespreid zijn. Er moet, met andere woorden, evenveel tijd tussen de observaties zijn. De temperatuurgegevens bevatten geen ontbrekende waarden, dus dit vormt geen probleem. Echter, de schrikkeljaren hebben steeds een dag meer dan de andere jaren, namelijk 29 februari. Wanneer we de observaties op de schrikkeldagen in de tijdreeks zouden houden, bekomen we een tijdreeks die niet regelmatig gespreid is. Sommige jaren tellen dan immers 366 dagen in plaats van 365. De makkelijkste oplossing om toch een tijdreeks te verkrijgen waarbij er evenveel tijd tussen de observaties zit, is om de observaties op de schrikkeldagen uit de dataset te verwijderen. Dit is een standaardmethode, we merken op dat Campbell and Diebold (2005), Alaton et al. (2002) en Benth and Benth (2007) dit ook deden om respectievelijk het model van Campbell en Diebold, het Alaton-model en het Benth-model te specifi¨eren. Wij hebben in het bijzonder de waarnemingen op 29 februari 1984, 1988, 1992, 1996, 2000, 2004 en 2008 uit de gegevens verwijderd. Op die manier telt elk jaar 365 dagen en kunnen we werken met 10 950 gegevens in verband met de dagtemperatuur. We zullen een model opbouwen gebruik makend van het statistisch programma R, de belangrijkste output wordt doorheen de tekst meegegeven, de gebruikte R-code vindt u in Bijlage B. De gegevens die ons door het KMI geleverd werden, stonden in een Excelfile. Om deze in R te kunnen inladen, hebben we het bestand opgeslagen als txt-file. Hierdoor konden we de gegevens zonder problemen in R inladen en aan het programma meegeven dat het deze gegevens als een tijdreeks moet behandelen. Om een model op te kunnen stellen, moeten we eerst een globaal beeld krijgen van de tijdreeks. Hiertoe bekijken we de plot in Figuur 2.3. In deze figuur is er duidelijk sprake van seizoenseffecten. Het is logisch dat we seizoenseffecten waarnemen in een tijdreeks die de dagtemperatuur beschrijft, in de zomer is het immers gemiddeld warmer dan in 8 Het Koninklijk Meteorologisch Instituut van Belgi¨e. 2.4. Praktijk 31 de winter. Een eventuele trend is niet onmiddellijk waarneembaar in Figuur 2.3. Om na Figuur 2.3: Plot van de tijdreeks. te gaan of er toch een (stijgende) trend aanwezig is, hebben we in Excel de gemiddelde jaartemperatuur9 in een grafiek uitgezet. Hierna hebben we aan de grafiek een lineaire trendlijn toegevoegd, zie Figuur 2.4. Deze lineaire trendlijn werd automatisch door Excel berekend, Excel maakt hiervoor gebruik van de methode van de kleinste kwadraten. Uit Figuur 2.4: Plot van de gemiddelde jaartemperatuur en een lineaire trendlijn. de figuur kunnen we opmaken dat er mogelijk een licht stijgende trend aanwezig is. We zouden voor deze trend ook een verklaring kunnen bedenken. We zouden de toename in 9 De gemiddelde jaartemperatuur werd berekend als het gemiddelde van de reeds berekende dagtem365 1 X peratuur over ´e´en jaar, namelijk via de formule Ti . 365 i=1 32 Hoofdstuk 2. Modellering van de temperatuur de gemiddelde jaartemperatuur eventueel kunnen toeschrijven aan de opwarming van de aarde, de verstedelijking,. . . Het is belangrijk om te weten of er seizoenseffecten en een trend aanwezig zijn in de tijdreeks, want de meeste modellen (AR, MA, ARMA) zijn gedefinieerd voor stationaire tijdreeksen (zie Definitie 6). Vooraleer we een model specifi¨eren, moeten we dus de trend en de seizoenseffecten uit de tijdreeks verwijderen. We zullen dit op twee manieren doen, door gebruik te maken van de differentie-operator enerzijds en door gebruik te maken van een Fourierreeks anderzijds. 2.4.4.1 Differentie-operator De differentie-operator ∆ wordt als volgt gedefinieerd ∆Yt = Yt − Yt−1 . Een lineaire trend kan worden uitgeschakeld door ∆ ´e´en keer toe te passen. Seizoenseffecten van orde s kunnen daarentegen ge¨elimineerd worden door het toepassen van de seizoensgebonden differentie-operator van orde s, ∆s Yt = Yt − Yt−s . In deze sectie zullen we dus de trend en seizoenseffecten op deze manier verwijderen uit de tijdreeks alvorens een AR-, MA- of ARMA-model te specifi¨eren. Het is altijd het eenvoudigst wanneer we een model kunnen specifi¨eren voor de oorspronkelijke tijdreeks, de resultaten zijn dan immers het makkelijkst te interpreteren. Wanneer we dit in het achterhoofd houden, is een SARIMA-model aangewezen. Het is een model dat op de oorspronkelijke tijdreeks toegepast wordt, maar waarbij er toch rekening wordt gehouden met de aanwezige trend en seizoenseffecten. Wanneer we via de differentie-operatoren de trend en seizoenseffecten verwijderen, zal er echter vaak nog enige seizoensgebondenheid in de correlatiestructuur van de tijdreeks achterblijven. Om te bepalen welk SARIMA-model het meest geschikt is voor de temperatuurtijdreeks, bekijken we het correlogram van de tijdreeks nadat we de differentie-operator en de seizoensgebonden differentie-operator op de tijdreeks hebben laten inwerken. Wanneer Yt de oorspronkelijke tijdreeks voorstelt, zien we in Figuur 2.5 het correlogram van de tijdreeks ∆∆365 Yt . Wanneer we inzoomen op het begin van het correlogram, krijgen we het resultaat in Figuur 2.6. Deze figuur suggereert een MA(7)-structuur, aangezien de eerste zeven autocorrelaties significant verschillend van nul zijn10 . In Figuur 2.5 merkten we eveneens significante correlaties op rond lag 365, dit is de orde van de seizoensgevoeligheid. Hierdoor is het vanzelfsprekend om een SARMA(0,7)(0,1)-model voor te stellen 10 Zie Sectie 2.1.1 voor meer informatie over de relatie tussen een MA-model en het correlogram. 2.4. Praktijk 33 Figuur 2.5: Correlogram van de tijdreeks ∆∆365 Yt . Figuur 2.6: Correlogram van de tijdreeks ∆∆365 Yt . voor de tijdreeks ∆∆365 Yt , dit komt overeen met een SARIMA(0,1,7)(0,1,1)-model voor de oorspronkelijke tijdreeks Yt . Wanneer we dit model echter proberen te specifi¨eren in R, stuiten we op een fout: maximum supported lag is 350. Deze fout geeft aan dat het niet zo gebruikelijk is om voor dagelijkse gegevens de seizoenseffecten te verwijderen via een differentie-operator. Door gebruik te maken van die differentie-operator gaan we een bepaalde dag vergelijken met dezelfde dag een jaar voordien. Maar waarom zou de temperatuur op 24/10/2013 enkel vergelijkbaar zijn met 24/10/2012 en niet met de temperatuur op elke andere dag in de maand oktober? Aangezien de SARIMA-modellen enkel deze vergelijking maken, worden ze zeer zelden gebruikt om dagelijkse tijdreeksen met seizoenseffecten te beschrijven. Een oplossing wordt gegeven door de seizoenseffec- 34 Hoofdstuk 2. Modellering van de temperatuur ten te modelleren via een Fourierreeks, zie Sectie 2.4.4.2. Wanneer we toch een eenvoudig, eerste model willen specifi¨eren, kunnen we om deze “fout” in R heen werken. We kunnen een AR-, MA- of ARMA-model proberen specifi¨eren voor de tijdreeks ∆∆365 Yt in plaats van rechtstreeks voor de oorspronkelijke tijdreeks. Hiertoe bekijken we ook eens het partieel correlogram van deze tijdreeks in Figuur 2.7, bijna alle zichtbare autocorrelaties zijn significant verschillend van nul. Figuur 2.7: Partieel correlogram van de tijdreeks ∆∆365 Yt . Via Figuur 2.6 is een MA(7)-model voor de tijdreeks ∆∆365 Yt een eerste voorstel. Wanneer we dit model schatten in R, het commando maakt gebruik van de maximumlikelihoodmethode om de parameters te schatten (we nemen normaliteit van de innovaties aan), krijgen we onderstaande output Series: temperatuurzondertrendenseizoen ARIMA(0,0,7) with non-zero mean Coefficients: s.e. ma1 ma2 ma3 ma4 ma5 ma6 ma7 0.0015 -0.2603 -0.2058 -0.1644 -0.1159 -0.1254 -0.0952 0.0097 0.0096 0.0096 0.0105 0.0102 0.0092 0.0094 intercept -2e-04 s.e. 9e-04 2.4. Praktijk 35 sigma^2 estimated as 6.904: AIC=50505.87 AICc=50505.89 log likelihood=-25243.94 BIC=50571.28 Om na te gaan of de termen in het model significant verschillend van nul zijn, controleren we of nul al dan niet een element is van het betrouwbaarheidsinterval [Zn − 1.96 · SE(Zn ); Zn + 1.96 · SE(Zn )], waarbij Zn een schatter en SE(Zn ) de standaardfout van die schatter is. Voor de ma1-term krijgen we bijvoorbeeld het betrouwbaarheidsinterval [0.0015 − 1.96 · 0.0097; 0.0015 + 1.96 · 0.0097] = [−0.017512; 0.020512], er geldt 0 ∈ [−0.017512; 0.020512] waaruit volgt dat de ma1-term niet significant is. Via dezelfde werkwijze vinden we dat alle andere termen wel significant zijn. Hierdoor houden we de ma1-term toch in het model. De AIC-waarde11 van dit model is 50 505.87, het is een in-sample criterium dat in tegenstelling tot de MSE ook penaliseert voor de complexiteit van het model wanneer men de goodness-of-fit van dat model test. Hierdoor is een kleine AIC-waarde beter. BIC en AICc testen eveneens voor de goodness-of-fit, er geldt ook dat een lage waarde voor deze criteria beter is. Nu rest ons nog dit model te valideren. Hiertoe gaan we na of de residuen van het gefitte model zich gedragen als witte ruis. We bekijken eerst en vooral het correlogram van de residuen. Daarnaast doen we ook de Box-Ljung-test aan de hand van de Q-statistiek kmax X ρˆ2k , deze statistiek wordt vaak gebruikt om te testen of de reeks afkomstig is Q=T k=1 van een witte ruis stochastisch proces. De nulhypothese van de test is dat de autocorrelaties van de residuen gezamenlijk gelijk aan nul zijn, H0 : ρ1 = ρ2 = . . . = ρkmax = 0. Over de keuze van kmax kan gediscussieerd worden, wij kiezen kmax = 20. Wanneer we deze test uitvoeren voor het gefitte MA(7)-model, verkrijgen we een p-waarde die kleiner is dan 2.2 · 10−16 . Er geldt p-waarde < 0.05, waardoor we de nulhypothese verwerpen op het significantieniveau 0.05. In het correlogram van de residuen, Figuur C.1 in Bijlage C, zien we ook dat enkele autocorrelaties significant verschillen van nul. Het MA(7)-model is met andere woorden niet geschikt om de tijdreeks te beschrijven. We gaan nu op zoek naar een ander model. Uit Figuur 2.7 kunnen we afleiden dat een gewoon AR-model geen goede keuze is12 . Hierdoor stellen we het eenvoudigste ARMAmodel voor, namelijk het ARMA(1,1)-model nog steeds voor de tijdreeks ∆∆365 Yt . Analoog aan het bovenstaande kunnen we het model gaan schatten, zowel de ar1- als de ma-1 term blijken significant. De AIC-waarde voor dit model is 51 511.28. Wanneer we testen of de residuen zich gedragen als witte ruis, krijgen we opnieuw een p-waarde die 11 Het Akaike-informatieccriterium ln(MSE) + 2 Tp met p het aantal parameters in het model, zie ook Sectie 2.4.1. 12 Zie Sectie 2.1.2 voor meer informatie over de relatie tussen een AR-model en het partieel correlogram. 36 Hoofdstuk 2. Modellering van de temperatuur kleiner is dan 2.2 · 10−16 . In het correlogram, zie Figuur C.2, zijn ook opnieuw enkele significante autocorrelaties waarneembaar. Het voorgestelde ARMA(1,1)-model is dus ook niet gevalideerd. We zouden ook eens een ARMA(1,2)-model kunnen proberen. Wanneer we dit model schatten zijn de drie termen (ar1, ma1 en ma2) significant, de AIC-waarde is 50 350.1. Het correlogram van de residuen van het gefitte ARMA(1,2)-model is veelbelovend, er vallen slechts enkele significante autocorrelaties waar te nemen (zie Figuur C.3). De test voor witte ruis levert echter een p-waarde van 0.009466 < 0.05, we verwerpen de nulhypothese. Het model is opnieuw niet gevalideerd. Vervolgens specifi¨eren we een ARMA(2,2)-model. De geschatte ar1-, ma1- en ma2parameters zijn significant, de ar2-parameter is echter niet significant. Het correlogram van de residuen ziet er opnieuw redelijk goed uit, zie Figuur C.4. De Box-Ljung-test daarentegen levert een p-waarde gelijk aan 0.01431 op, we verwerpen de nulhypothese, wat het model niet geschikt maakt om de tijdreeks te beschrijven. Omdat in het voorgaande model de ar2-parameter niet significant verschillend van nul was, gaan we nu verder met het verhogen van het aantal MA-termen in het ARMA-model met ´e´en AR-term. Het volgende model dat we specifi¨eren, is daardoor het ARMA(1,3)model. Bij schatting van het model blijkt de ma3-term niet significant te zijn. De AIC-waarde is 50 350.01. Alhoewel er in het correlogram van de residuen (zie Figuur C.5) weinig significante autocorrelaties zichtbaar zijn, kunnen we op basis van de BoxLjung-test het model nog steeds niet valideren. De test leverde immers een p-waarde van 0.01551 op en dit is kleiner dan het door ons gekozen significantieniveau 0.05. Omdat we nog steeds geen gevalideerd model gevonden hebben, gaan we door met het verhogen van het aantal MA-termen, ook al was de ma3-term van het vorige model niet significant. We schatten dus het ARMA(1,4)-model. De ar1-, ma1- en ma2-term zijn duidelijk significant, de ma3-term is opnieuw niet significant en de bijkomende ma4-term is net wel significant. De AIC-waarde voor dit model is 50 348.03. Het correlogram van de residuen is te vinden in Figuur C.6. De p-waarde van de test voor witte ruis is 0.04086, we verwerpen de nulhypothese nog net op het significantieniveau 0.05. We kunnen het model dus niet valideren. We stellen nu het ARMA(1,5)-model voor. Wanneer we dit model schatten, krijgen we de volgende output Series: temperatuurzondertrendenseizoen ARIMA(1,0,5) with non-zero mean 2.4. Praktijk 37 Coefficients: s.e. ar1 ma1 ma2 ma3 ma4 ma5 intercept 0.8097 -0.8138 -0.2544 0.0168 0.0184 0.0331 0 0.0149 0.0179 0.0125 0.0133 0.0132 0.0125 0 sigma^2 estimated as 6.797: AIC=50343.29 AICc=50343.3 log likelihood=-25163.64 BIC=50401.42 De ma3- en ma4-term in dit model is niet significant, de vier andere termen zijn echter wel significant verschillend van nul. Omdat de ma5-term wel nog significant is, laten we de ma3- en ma4-term ook nog tot het model behoren. De AIC-waarde van dit geschatte model is 50 343.29. Wanneer we overgaan tot de validatie van het model bekijken we het correlogram in Figuur 2.8. We nemen drie autocorrelaties waar die Figuur 2.8: Correlogram van de residuen van het gefitte ARMA(1,5)-model voor de tijdreeks ∆∆365 Yt . significant verschillend van nul zijn. Wanneer we de Box-Ljung-test uitvoeren, krijgen we het volgende resultaat Box-Ljung test data: mymodel6$res X-squared = 25.3615, df = 20, p-value = 0.188 Voor de p-waarde geldt 0.188 > 0.05, we verwerpen met andere woorden de nulhypothese niet. We kunnen dus besluiten dat de eerste twintig autocorrelaties gezamenlijk gelijk 38 Hoofdstuk 2. Modellering van de temperatuur aan nul zijn, dit betekent dat de residuen van het gefitte ARMA(1,5)-model zich gedragen als witte ruis. Het model is hiermee gevalideerd en geschikt bevonden om de tijdreekst ∆∆365 Yt te beschrijven. Tot slot hebben we ook nog gecontroleerd of we een beter model bekomen door nog enkele extra MA-termen toe te voegen. We hebben het ARMA(1,6)-, het ARMA(1,7)-, het ARMA(1,8)- en ARMA(1,9)-model ook nog geschat. Echter de extra MA-termen die er voor deze modellen bijkomen, zijn niet meer significant en de AIC-waarde zakt ook niet onder deze van het ARMA(1,5)-model. Door extra MA-parameters toe te voegen, kunnen we niet beter doen dan het ARMA(1,5)-model. We concluderen dus dat het ARMA(1,5)-model geschikt is om de tijdreeks ∆∆365 Yt te beschrijven. Voor de definitie van een algemeen ARMA-proces verwijzen we naar Sectie 2.1.3. Wanneer we de geschatte parameters van het ARMA(1,5)-model voor de tijdreeks ∆∆365 Yt invullen, krijgen we de volgende specificatie ∆∆365 Yt =0.8097∆∆365 Yt−1 + ut − 0.8138ut−1 − 0.2544ut−2 + 0.0168ut−3 + 0.0184ut−4 + 0.0331ut−5 . We kunnen het volgende resultaat gebruiken om deze gelijkheid om te vormen ∆∆365 Yt = ∆(Yt − Yt−365 ) = Yt − Yt−1 − Yt−365 + Yt−366 , we krijgen dus Yt − Yt−1 − Yt−365 + Yt−366 = 0.8097(Yt−1 − Yt−2 − Yt−366 + Yt−367 ) + ut − 0.8138ut−1 − 0.2544ut−2 + 0.0168ut−3 + 0.0184ut−4 + 0.0331ut−5 = 0.8097Yt−1 − 0.8097Yt−2 − 0.8097Yt−366 + 0.8097Yt−367 + ut − 0.8138ut−1 − 0.2544ut−2 + 0.0168ut−3 + 0.0184ut−4 + 0.0331ut−5 . Een model voor de Belgische temperatuur wordt dus gegeven door Yt =1.8097Yt−1 − 0.8097Yt−2 + Yt−365 − 1.8097Yt−366 + 0.8097Yt−367 + ut − 0.8138ut−1 − 0.2544ut−2 + 0.0168ut−3 + 0.0184ut−4 + 0.0331ut−5 , waarbij Yt staat voor de temperatuur op dag t. 2.4. Praktijk 2.4.4.2 39 Fourierreeks In het algemeen zijn de seizoensversies van de ARIMA-modellen ontworpen voor korte periodes van seizoenseffecten, met name 12 voor maandelijkse gegevens of 4 voor kwartaalgegevens. Het gebruik van de seizoensgebonden differentie-operator van zeer hoge orde heeft echter niet veel nut. Voor dagelijkse gegevens (periode 365) vergelijken we dan immers wat er vandaag gebeurd is met wat er precies een jaar geleden heeft plaatsgevonden, dit is voor een tijdreeks die de temperatuur beschrijft niet zinvol. Voor dagelijkse data is het beter, of alleszins meer standaard, om het seizoenspatroon te modelleren met behulp van een Fourierreeks. We gaan in deze sectie dan ook van startmet het schatten van de seizoenseffecten door 2πs 2πs steeds enkele sin q - en cos q -termen, met s = 1, . . . , 365 en q = 1, 2, 3, . . ., 365 365 toe te voegen aan een model met enkel een lineaire trend. We krijgen uiteindelijk een model voor de seizoenseffecten dat er als volgt uitziet 2πs 2πs 4πs 4πs fs = a + b · s + c1 sin + c2 cos + c3 sin + c4 cos + ··· . 365 365 365 365 Het is aan ons om te bepalen hoeveel sin- en cos-termen we moeten toevoegen om een geschikt seizoenspatroon te bekomen. We zullen het model voor de seizoenseffecten kiezen aan de hand van de determinatieco¨effici¨ent R2 en de AIC-waarde. We kiezen het model waarvoor deze waarden “in evenwicht” zijn, met name de R2 waarde zo hoog mogelijk en de AIC-waarde zo laag mogelijk. Wanneer we immers enkel naar een hoge waarde voor de determinatieco¨effici¨ent zouden streven, zou dit in het voordeel zijn van de complexere modellen (deze met meer parameters). Wanneer we het eerste model met enkel een lineaire trend gaan schatten (fs = a + b · s), krijgen we via de R-commando’s in Bijlage B.3 de volgende output Call: lm(formula = temperatuur ~ oneyear) Residuals: Min 1Q Median 3Q Max -20.8615 -4.5287 0.2171 4.7750 17.7506 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.4209962 0.1221703 68.93 <2e-16 *** 40 oneyear Hoofdstuk 2. Modellering van de temperatuur 0.0118918 0.0005786 20.55 <2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 6.379 on 10948 degrees of freedom Multiple R-squared: 0.03716,Adjusted R-squared: F-statistic: 422.5 on 1 and 10948 DF, 0.03707 p-value: < 2.2e-16 Hieruit kunnen we afleiden dat zowel het intercept als de oneyear-co¨effici¨ent significant verschillend van nul zijn, aangezien respectievelijk 0 ∈ / [8.421 − 1.96 · 0.122; 8.421 + 1.96 · 0.122] = [8.182; 8.660] en 0 ∈ / [0.012−1.96·0.001; 0.012+1.96·0.001] = [0.010; 0.014]. We zouden dit ook kunnen concluderen via de respectievelijke p-waarden, die beide kleiner zijn dan 0.05. Via de p-waarde bij de F-statistiek verwerpen we de nulhypothese dat alle regressieparameters, behalve deze van het intercept, gelijk zijn aan nul. In dit model is het logisch dat we deze nulhypothese verwerpen, aangezien we reeds gezien hebben dat de oneyear-co¨effici¨ent significant is. De waarde voor de determinatieco¨effici¨ent is 0.03716, wat betekent dat 3.716% van de variatie in de temperatuur verklaard wordt door het model. Wanneer we de AIC-waarde berekenen, bekomen we 71 659.55. Tot slot plotten we ook nog eens het patroon dat we in dit model beschrijven, we verwachten een rechte en dit is ook wat we te zien krijgen in Figuur 2.9. Figuur 2.9: Patroon beschreven door het model fs = a + b · s. We gaan nu over naar een volgend model door aan hetvorige een eerste sin-term toe te 2πs . Bij het schatten van dit voegen, hierdoor komen we aan fs = a + b · s + c1 sin 365 model, bekomen we de volgende output 2.4. Praktijk 41 Call: lm(formula = temperatuur ~ oneyear + sin(sp)) Residuals: Min 1Q Median 3Q Max -23.374 -4.430 0.245 4.665 17.119 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 11.3304685 0.1756989 64.488 < 2e-16 *** oneyear -0.0040069 0.0009032 -4.436 9.24e-06 *** sin(sp) -3.0385199 0.1345852 -22.577 < 2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 6.236 on 10947 degrees of freedom Multiple R-squared: 0.07999,Adjusted R-squared: F-statistic: 475.9 on 2 and 10947 DF, 0.07983 p-value: < 2.2e-16 Zowel het intercept, de oneyear-co¨effici¨ent als de sin-term zijn significant, want de bijhorende p-waarden zijn respectievelijk < 2 · 10−16 , 9.24 · 10−6 en < 2 · 10−16 en dit is allemaal kleiner dan 0.05. Via de p-waarde bij de F-statistiek (< 2.2 · 10−16 ) verwerpen we opnieuw de nulhypothese dat alle regressieparameters, behalve deze van het intercept, gelijk zijn aan nul. Gezamenlijk zijn de regressieparameters dus significant verschillend van nul. De waarde voor de determinatieco¨effici¨ent is 7.999%, wat groter is dan bij het voorgaande model. De AIC-waarde voor dit model is 71 163.21, dit is kleiner dan voor het voorgaande model. Op basis van de R2 - en AIC-waarde is dit dus een beter model. Ter vergelijking plotten we in Figuur 2.10 het patroon van dit model. Het voorgaande model was reeds beter dan het eerste, daardoor voegen wenu ook nog 2πs 2πs een extra cos-term toe en bekijken we fs = a + b · s + c1 sin + c2 cos . We 365 365 geven de verkregen output niet mee, maar bespreken deze wel kort. We zien dat zowel het intercept, de oneyear-co¨effici¨ent, de sin-term als de cos-term significant zijn, want de bijhorende p-waarden zijn respectievelijk < 2·10−16 , 1.88·10−10 , < 2·10−16 en < 2·10−16 . Via de p-waarde bij de F-statistiek (2.2 · 10−16 ) kunnen we opnieuw concluderen dat de regressieparameters gezamenlijk significant verschillend van nul zijn. De R2 -waarde is 72.31%, wat veel groter is dan bij het model met enkel een sin-term. De AIC-waarde 42 Hoofdstuk 2. Modellering van de temperatuur Figuur 2.10: Patroon beschreven door het model fs = a + b · s + c1 sin 2πs 365 . voor het model is 58 017.51, dit is op zijn beurt veel kleiner dan de AIC-waarde van het voorgaande model. We verkrijgen dus opnieuw een beter model. Vervolgens voegen weopnieuw toe, we bekomen fs = a + b · s + een extra sin-term 2πs 2πs 4πs c1 sin + c2 cos + c3 sin . Voor dit model is de oneyear-term niet 365 365 365 meer significant verschillend van nul, alle andere termen zijn wel significant. Echter via de p-waarde bij de F-statistiek (< 2.2 · 10−16 ) verwerpen we de nulhypothese dat alle regressieparameters, behalve deze van het intercept, gelijk zijn aan nul. We hoeven dus niet onmiddellijk de oneyear-term uit het model te verwijderen. De waarde voor de determinatieco¨effici¨ent is 72.4% en de AIC-waarde is 57 984.51. We merken op dat de R2 -waarde niet veel stijgt maar dat de AIC-waarde wel nog afneemt, we beschouwen dit dus als een beter model dan het voorgaande. We voegen toe enbekomen fs = a + b · aanhet voorgaande model nog eencos-term 2πs 2πs 4πs 4πs s + c1 sin + c2 cos + c3 sin + c4 cos . Voor dit model zijn 365 365 365 365 alle termen, behalve de oneyear-term, significant. De p-waarde bij de F-statistiek levert echter opnieuw een waarde die kleiner is dan 0.05, waardoor we de oneyear-term niet onmiddellijk laten vallen. De waarde van de determinatieco¨effici¨ent is in dit geval 72.42% en de AIC-waarde is 57 977.92. Omdat de AIC-waarde nog steeds daalt, ook al stijgt de R2 -waarde maar lichtjes, beschouwen we dit model opnieuw als een beter resultaat. We voegen opnieuween extra sin-term model toe om toteen volgend tekomen, namelijk 2πs 2πs 4πs 4πs 6πs +c2 cos +c3 sin +c4 cos +c5 sin . fs = a+b·s+c1 sin 365 365 365 365 365 Voor dit patroon zijn opnieuw alle termen significant, zowel afzonderlijk als gezamenlijk. De R2 -waarde is 72.44% en de AIC-waarde is 57 973.62, wat opnieuw respectievelijk iets 2.4. Praktijk 43 hoger en iets lager is dan voor het voorgaande model. We beschouwen dit model dus opnieuw als een betere voorstelling van de seizoenseffecten. Vervolgens voegen we een cos-term specificatie, we verkrijgen toe aande voorgaande 2πs 2πs 4πs 4πs dan fs = a + b · s + c1 sin + c2 cos + c3 sin + c4 cos + 365 365 365 365 6πs 6πs c5 sin + c6 cos . Voor dit model zijn alle termen, behalve de laatste cos365 365 term, significant verschillend van nul. De p-waarde bij de F-statistiek geeft echter aan dat de regressieparameters allen significant verschillend van nul zijn. Dit laatste is in het voordeel van onze niet significante cos-term. Maar we merken op dat de verkregen determinatieco¨effici¨ent gelijk is aan 72.44% en de AIC-waarde gelijk is aan 57 974.61. De R2 -waarde stijgt dus niet, terwijl de AIC-waarde wel groter wordt dan deze van het voorgaande model. Hieruit kunnen we besluiten dat dit model niet beter is dan het voorgaande. We kunnen op deze manier doorgaan en steeds een sin- en een cos-term toevoegen. Omdat de analyse van deze nieuwe modellen steeds analoog is, laten we dit hier buiten beschouwing. We geven onmiddellijk de output en analyse van het, volgens ons, meest geschikte model mee. Het beste model13 voldoet aan de volgende specificatie voor de seizoenseffecten fs = a + 8 X i=1 2πs 2πs 18πs cs,i sin i + cc,i cos i + cs,9 sin . 365 365 365 We merken op dat er geen lineaire term in s meer aanwezig is. We hebben deze term uit het model gelaten, omdat hij niet significant verschillend van nul was en het model zonder deze term een betere (lagere) AIC-waarde opleverde. De output voor deze specificatie is als volgt Call: lm(formula = temperatuur ~ sin(sp) + cos(sp) + sin(2 * sp) + cos(2 * sp) + sin(3 * sp) + cos(3 * sp) + sin(4 * sp) + cos(4 * sp) + sin(5 * sp) + cos(5 * sp) + sin(6 * sp) + cos(6 * sp) + sin(7 * sp) + cos(7 * sp) + sin(8 * sp) + cos(8 * sp) + sin(9 * sp)) 13 We merken hier op dat het geen verschil maakt of we eerst de sin- en dan pas de cos-term toevoegen 2πs of omgekeerd. We hebben immers steeds het koppel sin i 2πs 365 en cos i 365 , i = 1, . . . , 8, in het model opgenomen. Daarnaast hebben we ook getest of we een beter model konden bekomen door de sin 9 2πs 365 te vervangen door de corresponderende cos-term, dit leverde echter geen beter resultaat. Het koppel 2πs sin 9 2πs 365 en cos 9 365 leverde ook geen beter model op. 44 Hoofdstuk 2. Modellering van de temperatuur Residuals: Min 1Q Median 3Q Max -15.2086 -2.3178 -0.1242 2.3046 11.7484 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.597201 0.032460 326.473 < 2e-16 *** sin(sp) -2.572994 0.045905 -56.050 < 2e-16 *** cos(sp) -7.375602 0.045905 -160.671 < 2e-16 *** 0.397405 0.045905 8.657 < 2e-16 *** cos(2 * sp) -0.136062 0.045905 -2.964 0.00304 ** sin(3 * sp) -0.066995 0.045905 -1.459 0.14448 cos(3 * sp) 0.044629 0.045905 0.972 0.33097 sin(4 * sp) 0.250385 0.045905 5.454 5.02e-08 *** cos(4 * sp) 0.012405 0.045905 0.270 0.78698 sin(5 * sp) 0.029808 0.045905 0.649 0.51614 cos(5 * sp) 0.254952 0.045905 5.554 2.86e-08 *** sin(6 * sp) 0.184267 0.045905 4.014 6.01e-05 *** cos(6 * sp) -0.077406 0.045905 -1.686 0.09178 . sin(7 * sp) 0.078660 0.045905 1.714 0.08664 . cos(7 * sp) -0.185825 0.045905 -4.048 5.20e-05 *** sin(8 * sp) -0.206164 0.045905 -4.491 7.16e-06 *** cos(8 * sp) 0.006139 0.045905 0.134 sin(9 * sp) -0.126467 0.045905 -2.755 sin(2 * sp) 0.89361 0.00588 ** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 3.397 on 10932 degrees of freedom Multiple R-squared: F-statistic: 0.7274,Adjusted R-squared: 0.727 1716 on 17 and 10932 DF, p-value: < 2.2e-16 2πs 2πs 2πs 2πs We kunnen zien dat de sin 3 -, de cos 3 -, de cos 4 -, de sin 5 -, 365 365 365 365 2πs 2πs 2πs de cos 6 -, de sin 7 - en de cos 8 -term niet significant is. Doorheen 365 365 365 de zoektocht naar ons model hebben we toch besloten om deze termen in het model te 2.4. Praktijk 45 houden, omdat we ook betere AIC- en R2 -waarden verkregen zonder de termen te verwijderen. Daarnaast verwerpen we door de grootte van de p-waarde van de F-statistiek de nulhypothese die stelt dat alle regressieparameters, behalve deze van het intercept, gelijk zijn aan nul. De waarde voor de determinatieco¨effci¨ent is 0.7274, wat betekent dat 72.74% van de variatie in de temperatuur verklaard wordt door het model. Wanneer we de AIC-waarde berekenen, bekomen we 57 873.87. Tot slot plotten we in Figuur 2.11 het patroon dat via dit model beschreven wordt. Figuur 2.11: Patroon beschreven door het model fs 18πs cc,i cos i 2πs 365 ] + sin 365 . = a + P8 i=1 [cs,i sin i 2πs 365 + De R2 -waarde is hoog (72.74%), het model voor de seizoenseffecten is dus geschikt. Hierdoor zullen we onze analyse verder zetten met de residuen in plaats van met de geobserveerde tijdreeks voor de temperatuur. Voor deze residuen geldt residuen = temperatuur − seizoenseffecten (2.11) = temperatuur − fs . Verder werken met deze residuen is geen verlies van algemeenheid, want eens we het model kennen voor de residuen kunnen we het geschatte signaal fs er bij optellen om zo uiteindelijk het totale model voor de temperatuur te bekomen. Op deze residuen gaan we een gelijkaardige analyse toepassen als in Sectie 2.4.4.1, daar vindt u ook meer uitleg over de methodes die we hieronder toepassen. We gaan dus voor de residuen een AR-, MA- of ARMA-model proberen schatten. Om dit te doen, bekijken we in Figuur 2.12 eerst en vooral een plot van de tijdreeks die we restemp noemen, het is de tijdreeks die bestaat uit de residuen (2.11) waarvoor we een model proberen te schatten. 46 Hoofdstuk 2. Modellering van de temperatuur Figuur 2.12: Plot van de tijdreeks restemp. Uit deze figuur kunnen we afleiden dat dit een stationaire tijdreeks is. Dit is ook wel logisch, want we hebben de seizoenseffecten reeds verwijderd door enkel nog naar de residuen te kijken. Hierdoor kunnen we onmiddellijk verder gaan met onze analyse. We bekijken zowel het correlogram in Figuur 2.13 als het partieel correlogram in Figuur 2.14. In Figuur C.7 en Figuur C.8 in Bijlage C werd een uitgebreider correlogram en partieel correlogram opgenomen. Aan de hand van deze figuren, vooral het partieel correlogram, kunnen we een AR(3)-model voorstellen voor de tijdreeks restemp. Figuur 2.13: Correlogram van de tijdreeks restemp. 2.4. Praktijk 47 Figuur 2.14: Partieel correlogram van de tijdreeks restemp. We schatten dit AR(3)-model met behulp van R en krijgen de volgende output Series: restemp ARIMA(3,0,0) with non-zero mean Coefficients: s.e. ar1 ar2 ar3 intercept 1.0041 -0.2570 0.0672 0.0000 0.0095 0.0133 0.0095 0.0947 sigma^2 estimated as 3.387: AIC=44445.19 AICc=44445.19 log likelihood=-22217.59 BIC=44481.69 Via de constructie van het betrouwbaarheidsinterval van de ar1-term bekomen we 0 ∈ / [0.98548; 1.02272], waaruit blijkt dat de ar1-term significant verschillend van nul is. Analoog bekomen we voor de ar2-term en de ar3-term 0 ∈ / [−0.283068; −0.230932] en 0∈ / [0.04858; 0.08582], respectievelijk. De drie AR-termen zijn dus significant, het intercept is duidelijk niet significant verschillend van nul. Een verklaring hiervoor kan gevonden worden in de manier waarop we de seizoenseffecten hebben voorgesteld, daar was het intercept wel significant verschillend van nul. De AIC-waarde voor dit model is 44 445.19. Tot slot willen we dit model nog gevalideerd zien. Hiertoe gaan we na of de residuen van het gefitte AR(3)-model zich gedragen als witte ruis. We bekijken 48 Hoofdstuk 2. Modellering van de temperatuur om te beginnen het correlogram van deze residuen14 in Figuur 2.15. Op basis van dit Figuur 2.15: Correlogram van de residuen van het gefitte AR(3)-model voor de tijdreeks restemp. correlogram zouden we kunnen besluiten dat de residuen zich gedragen als witte ruis, er zijn immers slechts twee significante autocorrelaties waar te nemen. Om echt zeker te zijn dat de reeks afkomstig is van een witte ruis stochastisch proces, doen we ook nog de Box-Ljung-test. De output van deze test is als volgt Box-Ljung test data: mymodel$res X-squared = 22.5036, df = 20, p-value = 0.3138 De nulhypothese van deze test is dat de eerste twintig autocorrelaties van de residuen gezamenlijk gelijk zijn aan nul. Deze nulhypothese wordt hier duidelijk niet verworpen aangezien voor de p-waarde geldt 0.3138 > 0.05. We kunnen dus concluderen dat de residuen zich gedragen als witte ruis, het model is gevalideerd. Nu we een gevalideerd model gevonden hebben, zouden we het hierbij kunnen laten. Maar het zou kunnen dat het correlogram en het partieel correlogram van de tijdreeks restemp ons niet op weg gezet hebben naar het best mogelijke model, namelijk dat met de laagste AIC-waarde. Daarom stellen we ook nog enkele andere modellen voor, die we dan eveneens schatten en proberen valideren. 14 We willen nog eens duidelijk maken dat dit de residuen van het gefitte AR(3)-model zijn en niet de residuen waarmee we de tijdreeks restemp geconstrueerd hebben. 2.4. Praktijk 49 Voor de zekerheid hebben we ook eens geprobeerd om een AR(1)- en een AR(2)-model te schatten en te valideren. Maar zoals het partieel correlogram in Figuur 2.14 reeds deed vermoeden, kunnen we uit de Box-Ljung-test voor deze modellen besluiten dat dit geen geschikte modellen zijn voor de tijdreeks. We zouden aan het AR(3)-model nog een extra parameter kunnen toevoegen en het AR(4)-model overwegen. Uit de schatting van dit model blijkt dat de toegevoegde ar4-term niet significant verschillend van nul is. De AIC-waarde voor het AR(4)-model is 44 446.28, wat hoger is dan voor het AR(3)model. We besluiten dus dat het AR(4)-model niet beter is om de tijdreeks restemp te beschrijven. We hebben voor de zekerheid ook nog enkele AR-termen van hogere orde toegevoegd. Enkel het AR(5)-model doet iets beter dan het AR(3)-model, vanaf het AR(6)-model zijn de bijkomende termen niet meer significant. We bespreken hier nog kort de output van het AR(5)-model. Alle termen in dit model, behalve het intercept en de ar4-term, zijn significant. De AIC-waarde die verkregen wordt, is gelijk aan 44 443.76 en dit is iets lager dan de AIC-waarde van het AR(3)-model. We kunnen het AR(5)-model ook als geschikt beschouwen want de p-waarde van de Box-Ljung-test is 0.6446, waardoor we de nulhypothese niet verwerpen. Het correlogram van de residuen van het gefitte AR(5)-model is te vinden in Figuur C.9. Zoals het correlogram in Figuur 2.13 reeds impliceerde, zijn de MA-modellen ook niet geschikt om de tijdreeks restemp te beschrijven. We hebben dit ook nog eens gecontroleerd en kwamen tot de conclusie dat het eerste gevalideerde MA-model het MA(12)-model is. Dit model bevat reeds zeven parameters meer dan het AR(5)-model en levert ook (niet totaal onverwacht) een hogere AIC-waarde gelijk aan 44 463.94. Hierdoor verlaten we het pad van de MA-modellen. De enige overblijvende mogelijkheid is nu een beter model te vinden binnen de klasse van de ARMA-modellen. We gaan van start met het meest eenvoudige ARMA-model, namelijk het ARMA(1,1)-model. Wanneer we het ARMA(1,1)-model schatten, bekomen we in R de volgende output Series: restemp ARIMA(1,0,1) with non-zero mean Coefficients: s.e. ar1 ma1 intercept 0.7601 0.2462 0.0000 0.0073 0.0110 0.0914 50 Hoofdstuk 2. Modellering van de temperatuur sigma^2 estimated as 3.389: AIC=44450.4 log likelihood=-22221.2 AICc=44450.41 BIC=44479.61 Via de constructie van het betrouwbaarheidsinterval zien we dat zowel de ar1- als de ma1term significant zijn. De AIC-waarde voor het ARMA(1,1)-model is gelijk aan 44 450.4, wat hoger is dan deze van het AR(5)-model. De p-waarde die bij de Box-Ljung-test hoort, is zodanig dat het model net gevalideerd is, immers 0.07895 > 0.05. Het correlogram van het gefitte model is te vinden in Figuur C.10. Op basis van de AIC-waarde is het model niet beter dan het AR(5)-model, maar het doet het ook niet zo veel slechter. We moeten er immers rekening mee houden dat er bij de schatting van het AR(5)-model drie parameters extra moeten geschat worden. Toch geven we nog steeds de voorkeur aan het AR(5)-model, maar we bekijken ook nog enkele complexere ARMA-modellen, want misschien vinden we in die klasse van modellen nog een betere specificatie voor de tijdreeks restemp. Binnen de klasse van ARMA-modellen vinden we enkele modellen die op basis van de AIC-waarde beter presteren dan het AR(5)-model. Het beste model onder deze ARMAmodellen is het ARMA(2,2)-model. Wanneer we dit model schatten, krijgen we de volgende output Series: restemp ARIMA(2,0,2) with non-zero mean Coefficients: s.e. ar1 ar2 ma1 ma2 intercept 1.5584 -0.5971 -0.5569 -0.2189 -0.0008 0.1110 0.0869 0.1102 0.0243 0.1017 sigma^2 estimated as 3.386: AIC=44442.01 AICc=44442.01 log likelihood=-22215 BIC=44485.81 We zien dat alle termen, behalve het intercept, significant verschillen van nul via de constructie van het betrouwbaarheidsinterval [Zn − 1.96 · SE(Zn ); Zn + 1.96 · SE(Zn )], waarbij Zn een schatter en SE(Zn ) de standaardfout van die schatter is. We merken ook dat de AIC-waarde gelijk is aan 44 442.01 en dus lager dan deze van het AR(5)-model. Als het ARMA(2,2)-model nu ook nog gevalideerd is, beschouwen we dit model als een beter model. We bekijken in Figuur 2.16 het correlogram van de residuen. De residuen lijken zich te gedragen als witte ruis, er is slechts ´e´en significante autocorrelatie. 2.4. Praktijk 51 Figuur 2.16: Correlogram van de residuen van het gefitte ARMA(2,2)-model voor de tijdreeks restemp. Om helemaal zeker te zijn doen we de Box-Ljung-test, de output is als volgt Box-Ljung test data: mymodel6$res X-squared = 17.9023, df = 20, p-value = 0.5938 Voor de p-waarde geldt 0.5938 > 0.05, waardoor we de nulhypothese niet verwerpen. Het ARMA(2,2)-model is dus geschikt om de tijdreeks restemp te beschrijven. Omdat we in dit model een parameter minder moeten schatten dan in het AR(5)-model en de AICwaarde lager is, verkiezen we het ARMA(2,2)-model als specificatie voor de tijdreeks restemp. Nu we een model gevonden hebben voor de tijdreeks restemp, kunnen we het model voor de tijdreeks van de temperatuur opstellen. Dit model heeft de volgende vorm temperatuur = fourier seizoenseffecten + ARMA(2, 2). (2.12) Wanneer we vergelijking (2.12) invullen met wat we in het bovenstaande verkregen hebben, dan bekomen we het volgende model 8 X 2πt 2πt 18πt Yt =a + cs,i sin i + cc,i cos i + cs,9 sin 365 365 365 i=1 (2.13) + θ + α1 Yt−1 + α2 Yt−2 + ut + β1 ut−1 + β2 ut−2 , waar we Yt noteren voor de temperatuur op dag t. Om de notatie niet onnodig te verzwaren, zullen we de geschatte parameters niet in de vergelijking substitueren, maar 52 Hoofdstuk 2. Modellering van de temperatuur weergeven in Tabel 2.1. Aan de hand van dit model zullen we in Sectie 3.4.3 de prijs van een fictieve optie gaan bepalen. parameter geschatte waarde a 10.597201 cs,1 −2.572994 cc,1 −7.375602 cs,2 0.397405 cc,2 −0.136062 cs,3 −0.066995 cc,3 0.044629 cs,4 0.250385 cc,4 0.012405 cs,5 0.029808 cc,5 0.254952 cs,6 0.184267 cc,6 −0.077406 cs,7 0.078660 cc,7 −0.185825 cs,8 −0.206164 cc,8 0.006139 cs,9 −0.126467 θ −0.0008 α1 1.5584 α2 −0.5971 β1 −0.5569 β2 −0.2189 Tabel 2.1: De geschatte waarde van de parameters in model (2.13). 2.5 Conclusie In dit hoofdstuk hebben we besproken hoe we de temperatuur kunnen modelleren. Met het oog op het prijzen van temperatuurderivaten worden er doorgaans twee groepen van 2.5. Conclusie 53 modellering beschouwd, de discrete modellering enerzijds en de continue modellering anderzijds. Met discrete modellering worden de modellen die de temperatuur beschrijven via econometrische processen bedoeld en met continue modellering worden de modellen die de temperatuur voorstellen via een stochastische differentiaalvergelijkingen bedoeld. In dit hoofdstuk zijn we dan ook begonnen met het geven van enkele standaard econometrische modellen en stochastische differentiaalvergelijkingen die in de literatuur vaak gebruikt worden om de temperatuur te beschrijven. Dit zijn bijvoorbeeld het AR-, MA-, ARMA-, Ornstein-Uhlenbeck- en Hull-White-proces. Vervolgens hebben we een voorbeeld gegeven van een discreet model dat in de literatuur vaak geciteerd wordt, dit is het model van Campbell en Diebold (zie Sectie 2.4.1). We hebben ook twee voorbeelden van een continu model voor de temperatuur gegeven, namelijk het Alaton-model en het Benth-model. Tot slot hebben we ook zelf twee modellen opgesteld voor de Belgische temperatuur. Deze twee modellen beschrijven de temperatuur op dag t, deze variabele wordt in de modellen voorgesteld door Yt . Het eerste model maakt gebruik van de differentie-operator om de seizoenseffecten te modelleren, dit resulteert in Yt =1.8097Yt−1 − 0.8097Yt−2 + Yt−365 − 1.8097Yt−366 + 0.8097Yt−367 + ut − 0.8138ut−1 − 0.2544ut−2 + 0.0168ut−3 + 0.0184ut−4 + 0.0331ut−5 . Het tweede model maakt daarentegen gebruik van een Fourierreeks om de seizoenseffecten te beschrijven, dit levert 8 X 2πt 2πt 18πt Yt =a + cs,i sin i + cc,i cos i + cs,9 sin 365 365 365 i=1 + θ + α1 Yt−1 + α2 Yt−2 + ut + β1 ut−1 + β2 ut−2 , waarbij de numerieke waarde voor de parameters opgelijst staat in Tabel 2.1. 54 Hoofdstuk 2. Modellering van de temperatuur Hoofdstuk 3 Temperatuurderivaten prijzen In dit hoofdstuk zullen we het prijzen van temperatuurderivaten bespreken. De prijs van zo een derivaat zal afhankelijk zijn van het type contract, de temperatuurindex,. . . Hier vinden we dan ook de link met Hoofdstuk 2. Aangezien de waarde van het derivaat afhankelijk is van de gekozen temperatuurindex, is het noodzakelijk dat we de temperatuur kunnen voorspellen. Hiertoe gebruiken we de modellen gedefinieerd in vorige sectie. Het prijzen van het derivaat zal dus ook sterk afhangen van het gebruikte model en de correctheid ervan. Wanneer men aandelenopties wil gaan prijzen, doet men beroep op het Black-ScholesMerton-model, we zullen zien dat voor temperatuurderivaten zo een algemeen model niet bestaat. Er bestaat dus nog geen algemeen aanvaard, zowel praktisch bruikbaar als wiskundig correct, kader voor het prijzen van temperatuurderivaten (Alexandridis and Zapranis (2013)). Echter, in de praktijk merken we dat de temperatuurderivaten steeds vaker verhandeld worden, en er dus toch enkele gangbare methoden moeten bestaan om ze te prijzen. In dit hoofdstuk zullen we deze methoden en hun wiskundige correctheid belichten. We gaan van start met het bespreken van het bekende Black-Scholes-Mertonmodel en in het bijzonder de tekortkomingen ervan wanneer we het willen gebruiken om temperatuurderivaten te prijzen. Vervolgens bespreken we drie methoden, om temperatuurderivaten te prijzen, die in de literatuur het meest behandeld worden. De eenvoudigste methode is historical burn analysis. We bespreken deze methode en berekenen de prijs van een fictieve calloptie voor Belgische temperatuurgegevens in Sectie 3.2. In Sectie 3.3 bespreken we een iets nauwkeurigere methode, namelijk indexmodellering. De meest nauwkeurige methode is dagelijkse simulatie. Deze methode maakt gebruik van een model voor de temperatuur. Op basis van ´e´en van de econometrische modellen voor de Belgische temperatuur die we in Hoofdstuk 2 hebben opgesteld, tonen we hoe een 55 56 Hoofdstuk 3. Temperatuurderivaten prijzen optie geprijsd kan worden via discrete dagelijkse simulatie. In Sectie 3.4.4 tonen we ook aan hoe men via continue dagelijkse simulatie te werk gaat om specifieke temperatuurderivaten (futures en callopties) te prijzen. Tot slot bespreken we ook nog twee eerder alternatieve methoden om temperatuurderivaten te prijzen, beide methoden maken gebruik van zogenaamde nutsfuncties. 3.1 Black-Scholes-Merton Aangezien we op zoek zijn naar een methode voor het prijzen van temperatuurderivaten, is het een logische eerste stap om eens te kijken naar het beroemde Black-Scholes-Mertonmodel. Dit model vormt immers de basis voor het prijzen van andere financi¨ele derivaten. Zoals in Shreve (2004) kan de Black-Scholes-Merton-differentiaalvergelijking opgesteld worden, dit is een parti¨ele differentiaalvergelijking voor de prijs van een optie waarbij het prijsproces van het onderliggend aandeel gemodelleerd is als een geometrische Brownse beweging. De gezochte continue functie c(t, x) voor de waarde van de calloptie is immers de oplossing van de Black-Scholes-Merton parti¨ele differentiaalvergelijking 1 ct (t, x) + rxcx (t, x) + σ 2 x2 cxx (t, x) = rc(t, x) voor elke t ∈ [0, T ), x ≥ 0, 2 (3.1) met eindvoorwaarde c(T, x) = (x − K)+ , waarbij (x − K)+ = max{x − K, 0}. Wanneer we naast de eindvoorwaarde ook enkele randvoorwaarden specificeren kan de parti¨ele differentiaalvergelijking opgelost worden. De oplossing is c(t, x) = xΦ(d+ (T − t, x)) − Ke−r(T −t) Φ(d− (T − t, x)), met 0 ≤ t < T, x > 0, 1 σ2 x d± (τ, x) = √ ln + r± τ , K 2 σ τ en Φ de standaard normale cumulatieve verdelingsfunctie, d.w.z. Z y z2 1 e− 2 dz. Φ(y) = √ 2π −∞ In de literatuur wordt er echter niet veel aandacht geschonken aan het gebruik van het Black-Scholes-Merton-model voor het prijzen van weerderivaten. Het Black-ScholesMerton-model toont immers wel hoe de prijs van de optie bepaald kan worden uit het onderliggend aandeel, maar maakt hiertoe ook enkele belangrijke onderstellingen: • Het aandeel kan continu verhandeld worden, zonder extra transactiekosten. 3.1. Black-Scholes-Merton 57 • De prijs van het aandeel volgt een geometrische Brownse beweging met drift en wordt niet be¨ınvloed door het verhandelen van het aandeel (deze transacties worden ondernomen om de optie te hedgen). • De markt werkt op zo een manier dat elke mogelijkheid tot arbitrage onmiddellijk verdwijnt. Bij het opstellen van de Black-Scholes-Merton parti¨ele differentiaalvergelijking gaat men er van uit dat de waarde van de call overeenkomt met de portefeuillewaarde (de waarde van een hedgingstrategie die belegt in de spaarrekening en in het aandeel) en er eveneens geen mogelijkheid tot arbitrage bestaat. Het Black-Scholes-Merton-model is een houvast bij het prijzen van derivaten in een complete markt. Een financi¨ele markt is compleet wanneer elk afgeleid product kan gehedged worden (Shreve (2004)). Desondanks het succes, brokkelt de aanpak van Black, Scholes en Merton af in incomplete financi¨ele markten, er is dan niet meer voldaan aan de drie bovenstaande onderstellingen (hoofdzakelijk de eerste). De markt voor weerderivaten is een voorbeeld van een incomplete markt1 , wat betekent dat niet elk product kan gehedged worden. Het grootste verschil tussen traditionele financi¨ele derivaten en weerderivaten is dat bij de laatstgenoemde de prijzen zijn gekoppeld aan het weer (bij temperatuurderivaten in het bijzonder aan de temperatuur) in plaats van aan een onderliggend aandeel. Zoals in Hoofdstuk 1 reeds vermeld werd, betekent dit dat de onderliggende geen waarde heeft en dat er dus gebruik gemaakt wordt van temperatuurindices als onderliggende. Deze indices (bijvoorbeeld HDD en CDD2 ) zijn niet verhandelbaar. Anderzijds is er ook nog weinig of geen liquiditeit in weerderivaten. Op die manier zijn de traditionele niet-arbitrage modellen voor het prijzen van financi¨ele derivaten, zoals het Black-Scholes-Merton-model, niet toepasbaar om weerderivaten te prijzen. Ten slotte vermelden we nog dat de papers die wel aandacht besteden aan het gebruik van het Black-Scholes-Merton-model voor het prijzen van temperatuurderivaten (zie bijvoorbeeld Meissner and Burke (2011) en Jewson and Zervos (2003)) eerder de aangepaste versie voor forwards of futures bespreken. Deze derivaten verschillen van opties in het feit dat opties aan de houder het recht geven om te kopen of te verkopen, terwijl futures en forwards de verplichting om te kopen of verkopen meebrengen. We merken tevens op dat wanneer de intrestvoet niet stochastisch is, de forward- en de futureprijs samenvallen 1 De markt voor het weer is een incomplete markt aangezien de onderliggende van het weerderivaat niet kan opgeslagen of verhandeld worden. 2 Zie Definitie 2. 58 Hoofdstuk 3. Temperatuurderivaten prijzen (cf. Shreve (2004)). Er wordt in de voornoemde papers gesteld dat wanneer er geen rekening gehouden wordt met het feit dat de financi¨ele markt voor temperatuurderivaten incompleet is, we gewoon gebruik kunnen maken van het Black-Scholes-Merton-model. Immers, de payoff van de opties (call of put) met als onderliggende een temperatuurindex kan op dezelfde manier gedefinieerd worden als deze van opties op aandelen. Het enige wat we dan echter in ons achterhoofd moeten houden, is dat de prijs van de future of de forward in dit bijzonder geval afhankelijk is van de temperatuurindex. Wat we moeten onthouden over deze uiteenzetting is dat het Black-Scholes-Mertonmodel in de praktijk niet gebruikt wordt om temperatuurderivaten te prijzen, omdat de markt voor weerderivaten incompleet is. 3.2 Historical burn analysis Het Black-Scholes-Merton-model wordt niet gebruikt in de context van temperatuurderivaten, hierdoor zijn er vele andere methoden ontstaan om temperatuurderivaten te prijzen. De eenvoudigste van deze methoden is Historical Burn Analysis3 (HBA). In de literatuur (zie bijvoorbeeld Alexandridis and Zapranis (2013) en Harris (2003)) wordt deze methode als volgt beschreven. HBA is een klassieke aanpak voor het prijzen van weerderivaten waarbij simulaties gebaseerd op historische data uitgevoerd worden. Meer bepaald, HBA berekent de gemiddelde payoff van de weerderivaten in de afgelopen n jaar. Deze methode waardeert de weerderivaten dus gebaseerd op de payoff die men zou verkregen hebben, wanneer het contract in het verleden werd gehouden. Men stelt zich dus de vraag “Wat zou de uitbetaling geweest zijn wanneer we een gelijkaardige optie de voorbije n jaar elk jaar zouden hebben verkocht?”Vaak wordt n tussen 10 en 30 jaar gekozen. De keuze van n vereist een zorgvuldige beoordeling van het evenwicht tussen het gebruik van een tijdreeks die lang genoeg is om een redelijk niveau van statistische significantie te verkrijgen en een reeks die kort genoeg is om conclusies te kunnen trekken die nog voldoende relevant zijn. Wanneer we historical burn analysis willen toepassen in de praktijk, hanteren we de volgende werkwijze. Eerst en vooral moeten we historische temperatuurdata verzamelen en de noodzakelijke correcties aanbrengen. De correcties zijn nodig omdat de verkregen gegevens vaak enkele problemen vertonen: • Ook al is de data beschikbaar, deze kan foute gegevens bevatten en/of er kunnen zelfs gegevens ontbreken. 3 Deze methode staat ook bekend als actuarieel prijzen. 3.2. Historical burn analysis 59 • Een schrikkeljaar heeft 366 dagen in plaats van 365. • Het weerstation dat de gegevens levert, kan gedurende de meetperiode significant gewijzigd zijn. Het weerkundig observatorium kan doorheen de jaren bijvoorbeeld elders gevestigd worden of er kan bijvoorbeeld een stadseffect ontstaan (de bebouwing rondom het weerstation neemt toe en beton en bakstenen houden de warmte veel langer vast). • De gegevens zouden een bepaald jaar kunnen bevatten waarin extreme weersomstandigheden voorkwamen. Het is echter mogelijk om de data aan te passen zodat deze problemen omzeild worden. We kunnen dan eveneens de data omzetten naar de gebruikte index, zie Hoofdstuk 1. Vervolgens moeten we voor de voorbije n jaren bepalen wat de optie zou hebben uitbetaald en het gemiddelde van deze n payoffs berekenen. Ten slotte moeten we dit gemiddelde verdisconteren naar de waarderingsdatum (de dag waarop het contract begint te lopen). De prijs wordt dus n 1X P (t) = D(t, T ) · payoffi , n i=1 (3.2) waarbij D(t, T ) de verdisconteringsfactor van maturiteit T naar t voorstelt. Er bestaat echter ook nog een onvoorspelbare component van de weerschommelingen, deze zit vervat in het weerrisico. Wanneer we dus naast het gemiddelde ook de standaarddeviatie van de payoffs beschouwen en het weerrisico voorstellen door deze standaarddeviatie, dan kunnen we de prijs van het contract als volgt neerschrijven P (t) = D(t, T ) · (µ ± α · σ), (3.3) waarbij D(t, T ) de verdisconteringsfactor van maturiteit T naar t voorstelt, µ de historische gemiddelde payoff is, σ de historische standaarddeviatie weergeeft en α een positief getal is dat de risicotolerantie uitdrukt. Het is duidelijk dat HBA zeer simplistisch is en ook zeer makkelijk te berekenen, het is immers niet nodig om de verdeling van de temperatuur te beschouwen of een stochastische differentiaalvergelijking op te lossen. Bovendien worden er in deze methode erg weinig onderstellingen gemaakt. De eerste en voornaamste onderstelling is dat de tijdreeks die de temperatuur weergeeft stationair is. Een tijdreeks is, eenvoudig gezegd, stationair (zie Definitie 6) wanneer ze voor elk tijdsvenster dezelfde algemene eigenschappen heeft. Dit houdt in dat de toekomst zich ongeveer zoals het verleden zal gedragen. 60 Hoofdstuk 3. Temperatuurderivaten prijzen Daarnaast wordt er ook nog aangenomen dat de data voor de verschillende jaren onafhankelijk en identiek verdeeld zijn. Door de eenvoud van deze methode zijn er ook verschillende nadelen aan verbonden. Eerst en vooral, wanneer we een tijdreeks van de temperatuur inspecteren, merken we dat de gemaakte onderstellingen eigenlijk niet correct zijn. Zo geldt stationariteit niet voor een tijdreeks die de temperatuur weergeeft. Studies tonen immers aan dat de temperatuur een stijgende trend vertoont (als gevolg van verstedelijking, opwarming van de aarde,. . . ) en er eveneens seizoensgevoeligheid aanwezig is. Dit hebben we in Sectie 2.4.4 ook opgemerkt, we verwijzen hiervoor in het bijzonder naar Figuur 2.3 en 2.4. Daarnaast geldt er voor observaties in een tijdreeks dat deze niet onafhankelijk zijn. Vervolgens zal de aanwezigheid van een afwijking in de temperatuurdata4 een significant effect hebben op de waardering van het weerderivaat. Omdat de onderstellingen niet correct zijn, er geen voorspellingen worden gebruikt en het resultaat sterk be¨ınvloed wordt door de gebruikte data, zal HBA een vertekende en onnauwkeurige methode zijn. Schiller et al. (2012) breiden het concept van de HBA-methode, zoals deze hierboven beschreven werd, uit. Ze stellen eveneens dat temperatuurderivaten gewaardeerd kunnen worden door de toekomstige payoff van het contract te bepalen uit de gegeven historische payoffs. Daarnaast benadrukken ze dat de index gekend dient te zijn om de payoff van het contract te bepalen. Wanneer men bijvoorbeeld een derivaat voor de meetperiode [τ1 , τ2 ] wil prijzen voor het jaar n + 1, kan men de (fictieve) indices (voor hetzelfde derivaat) in de jaren n, n − 1, n − 2,. . . gaan berekenen. Deze indices kunnen bepaald worden uit de gebruikte temperatuurdata. Hierdoor bekomt men een stochastisch proces dat de indices van de voorbije n jaren voorstelt: Y1 , Y2 ,. . . ,Yn . Wat Schiller et al. (2012) nu extra doen is op basis van deze informatie de index voor het jaar n + 1 (d.i. Yn+1 ) proberen voorspellen. Men kan als volgt te werk gaan. Vertrek van het simpel lineair regressiemodel Yi = β0 + β1 · i + εi , i = 1, . . . , n, waarbij • E[εi ] = 0, i = 1, . . . , n; • Var(εi ) = σ 2 , i = 1, . . . , n; • Cov(εi , εj ) = 0, i 6= j. 4 Deze afwijking kan veroorzaakt worden door zowel extreme weersomstandigheden als fouten in de data. 3.2. Historical burn analysis 61 Hieruit kan men β0 (intercept) en β1 (rico) via de methode van de kleinste kwadraten5 als volgt schatten (cf. Goetghebeur (2011)) ! n n X X 1 n+1 ˆ βˆ0 = · Yi − βˆ1 i =Y − · β1 , n 2 i=1 i=1 Pn Pn Pn n+1 i=1 i (i − )(Y − Y ) i i=1 i=1 (i − 2 )(Yi − Y ) n ˆ P P β1 = , = n Pn n n+1 2 i=1 i 2 ) i=1 (i − 2 ) i=1 (i − n n n X 1X n · (n + 1) Yn en i= . n i=1 2 i=1 Uit de Gauss-Markov stelling weten we dat onder de voorwaarden van het simpel lineair regressiemodel met ongekende foutterm-verdeling, de kleinste kwadraten schatters βˆ0 en met Y = βˆ1 onvertekend zijn en minimale variantie6 hebben. Samengevat betekent dit dat onder de voorwaarden E[εi ] = 0, Var(εi ) = σ 2 en Cov(εi , εj ) = 0, i 6= j, de beste lineaire onvertekende schatter voor Yi gelijk is aan Ybi = βˆ0 + βˆ1 · i. Dus, de index Yn+1 kan als volgt voorspeld worden Ybn+1 = βˆ0 + βˆ1 · (n + 1). Nu kan men dan de prijs van het derivaat opstellen op basis van de payoff die door gebruik van deze index bekomen wordt. Samengevat kan men stellen dat Schiller et al. (2012) niet onmiddellijk de historische payoffs beschouwen, maar wel de historische index. Op basis van de gegevens over de historische index gaan zij een model opstellen en gebruiken dit model om de toekomstige index te bepalen. Die index gebruiken ze dan om het derivaat te prijzen. We merken nog op dat de in de literatuur meest voorkomende beschrijving van HBA echter geen voorspellingen incorporeert. Historical burn analysis wordt beschouwd als de eenvoudigste methode (op het gebied van implementatie) om een derivaat te prijzen en als de meest gevoelige voor grote fouten in de berekende prijs. Toch wordt HBA nog steeds beschouwd als een aanvaardbare eerste benadering om de prijs van het contract te bepalen, hierdoor wordt het ook veel gebruikt door marktdeelnemers. Omdat deze methode zo eenvoudig is, gebruiken we ze om de prijs te bepalen van een calloptie met behulp van Belgische temperatuurdata. 5 De methode van de kleinste kwadraten vertrekt van het simpel lineair regressiemodel Yi = β0 + Pn 2 β1 · xi + εi en zoekt dewaarden voor β0 en β1 die i=1 (yi − β0 − β 1 xi ) minimaliseren. Deze worden Pn P P (x −¯ x)(Yi −Y ) n n ˆ ˆ ¯ en βˆ1 = i=1 Pn i gegeven door βˆ0 = n1 · . i=1 xi = Y − β1 x i=1 Yi − β1 x) 2 i=1 (xi −¯ P n 6 Elke onvertekende schatter van de vorm i=1 wi Yi heeft een variantie die minstens gelijk is aan de variantie van de kleinste kwadraten schatter. 62 Hoofdstuk 3. Temperatuurderivaten prijzen 3.2.1 Prijzen op basis van HBA met Belgische temperatuurdata In deze sectie gaan we via HBA een fictief contract7 prijzen gebruik makend van de Belgische temperatuurdata die we ook in Sectie 2.4.4 gebruikt hebben. Het weerstation waarvan de gegevens worden verkregen, is gevestigd in Ukkel. Het contract is een calloptie geschreven op de CAT-index voor de zomermaanden (in dit geval mei-september). Het contract heeft dus een looptijd van 153 dagen. We veronderstellen dat het uitoefenniveau 2 500 is, dat de tick size e 20 per indexpunt is en dat er geen maximum op de uitbetaling staat. We merken op dat de gebruikte temperatuurdata slechts lopen tot en met 2010 en daarom nemen we aan dat we dit contract voor de zomermaanden in 2011 willen kopen. Om de prijs van de calloptie via HBA te bepalen, gaan we eerst en vooral na wat de optie zou hebben uitbetaald in de voorbije jaren. De payoff van deze calloptie ziet er als volgt uit payoff = max{CAT − 2 500, 0} × tick size. (3.4) De berekeningen zijn eenvoudig in Excel te maken. We bepalen eerst voor elk jaar de CAT-index8 voor de periode mei-september zoals in Definitie 3. Uit de waarde van de index kunnen we vervolgens max{CAT − 2 500, 0} bepalen. Uiteindelijk kunnen we de payoff berekenen door deze waarde te vermenigvuldigen met e 20, zie Tabel 3.1. Wanneer de prijs van de calloptie bepaald wordt via (3.2), moeten we het gemiddelde van de payoffs berekenen. Een simpele calculatie levert dat dit gemiddelde gelijk is aan e 854. De prijs voor de calloptie op tijdstip t wordt dus P (t) = D(t, T ) · 854, (3.5) waarbij D(t, T ) de verdisconteringsfactor van maturiteit T naar t voorstelt. Wanneer we ook het weerrisico in onze prijs willen opnemen, dan kunnen we de prijs van de calloptie bepalen via (3.3). Hier hebben we ook de historische standaarddeviatie van de payoff voor nodig, deze waarde kan ook via een ingebouwde formule in Excel berekend worden. Deze calculatie levert σ = 1 275.69. De prijs voor de calloptie op tijdstip t wordt in dit geval P (t) = D(t, T ) · (854 ± α · 1 275.69), 7 8 De parameters van zo een contract zijn te vinden in Hoofdstuk 1. We ronden af op de eenheid, omdat de tick size gedefinieerd is per indexpunt. (3.6) 3.3. Indexmodellering 63 waarbij D(t, T ) de verdisconteringsfactor van maturiteit T naar t voorstelt en α een positief getal is dat de risicotolerantie uitdrukt. Omdat de prijs van de calloptie positief 854 moet zijn, moet 854 ± α · 1 275.69 ≥ 0 gelden, d.i. α ≤ = 0.67. 1 275.69 Wanneer we veronderstellen dat r = 2.25% de constante intrestvoet is en D(t, T ) = e−r(T −t) noteren, kunnen we via (3.2) de prijs voor de calloptie op de ingangsdatum van het contract bepalen. Formule (3.5) verder invullen, levert 153 P (t) = e−0.0225· 365 · 854 = 845.98. De prijs voor de calloptie op de ingangsdatum van het contract is dus e 845.98. Wanneer we nu ook nog veronderstellen dat voor de risicotolerantie α = 0.5 geldt, kunnen we via (3.6) grenzen voor de prijs van de calloptie bepalen. We krijgen 153 P (t) = e−0.0225· 365 · (854 ± 0.5 · 1 275.69) = e−0.0094 · (854 ± 637.845), en dus P (t) ∈ [214.13, 1 477.84]. De prijs voor de calloptie op de ingangsdatum van het contract varieert dus tussen e 214.13 en e 1477.84 wanneer we ook het weerrisico meenemen in de prijsbepaling. 3.3 Indexmodellering In deze sectie zullen we een methode bespreken die nog steeds redelijk eenvoudig is, maar vaak nauwkeurigere resultaten zal opleveren dan de aanpak via historical burn analysis, die in de vorige sectie uiteengezet werd. We baseren ons op Alexandridis and Zapranis (2013) en Brix et al. (2002). De methode vereist het opstellen van een model voor de onderliggende temperatuurindex (zie Hoofdstuk 1) en staat bekend als indexmodellering. Men kan de index modelleren op basis van zowel parametrische als niet-parametrische verdelingen. Een eerste stap in de modellering is de specificatie van het model, hiertoe moet men een verdeling kiezen. Brix et al. (2002) stellen dat de eenvoudigste vorm van een nietparametrische verdeling, de empirische (historische) distributie van de indices is. Het gebruik van deze distributie wordt ook burn analysis genoemd. Met andere woorden, modellering van de index via een niet-parametrische verdeling valt samen met de methode uit Sectie 3.2. Meer bepaald met de methode die door Schiller et al. (2012) beschreven werd, omdat deze specifiek focust op de temperatuurindices. Voor het verschil tussen HBA en indexmodellering volgens Schiller et al. (2012), zie infra. 64 Hoofdstuk 3. Temperatuurderivaten prijzen jaar CAT-index max{CAT − 2 500, 0} payoff (e) 1981 2 381 0 0 1982 2 583 83 1 660 1983 2 553 53 1 060 1984 2 282 0 0 1985 2 365 0 0 1986 2 348 0 0 1987 2 305 0 0 1988 2 390 0 0 1989 2 599 99 1 980 1990 2 489 0 0 1991 2 389 0 0 1992 2 568 68 1 360 1993 2 385 0 0 1994 2 527 27 540 1995 2 586 86 1 720 1996 2 295 0 0 1997 2 551 51 1 020 1998 2 474 0 0 1999 2 619 119 2 380 2000 2 498 0 0 2001 2 486 0 0 2002 2 512 12 240 2003 2 692 192 3 840 2004 2 499 0 0 2005 2 569 69 1 380 2006 2 756 256 5 120 2007 2 498 0 0 2008 2 547 47 940 2009 2 619 119 2 380 2010 2 497 0 0 Tabel 3.1: Gegegevens die de payoff van het contract bepalen. 3.3. Indexmodellering 65 Wanneer men het heeft over indexmodellering wordt dus eigenlijk modellering van de index via parametrische verdelingen bedoeld. Men kan verschillende parametrische verdelingen gebruiken en dit afhankelijk van de gebruikte temperatuurindex. Een normale verdeling is vaak een gepaste keuze voor de HDD en CDD indices, vooral voor contracten over een langere periode. Dit kan verklaard worden door de centrale limietstelling. Stelling 1. Als X1 , X2 , . . . , Xn onafhankelijk en identiek verdeelde stochastische variabelen zijn met gemiddelde µ en variantie 0 < σ 2 < ∞, dan geldt " n # X √ lim Pr Xi ≤ nµ + xσ n = Φ(x), n→∞ i=1 met Φ(x) de standaard normale verdelingsfunctie. Er kan ook worden aangetoond dat wanneer extreme gebeurtenissen (zoals het aantal dagen met een temperatuur die een hoge drempel overschrijdt) bijna onafhankelijk zijn, ze ongeveer Poisson verdeeld moeten zijn. Echter, extreme weersomstandigheden komen vaak voor in clusters en kunnen dus vaak beter gemodelleerd worden door het gebruik van een negatief binomiale verdeling. Er bestaat echter praktisch geen theorie die aangeeft welke verdeling effectief het best gebruikt wordt om de index te beschrijven. Hierdoor bestaat er een significante kans dat een niet geschikte verdeling gebruikt wordt, waardoor er grote fouten ontstaan in de verkregen prijzen van het contract. Een volgende stap in de modellering is het schatten van het model. In de parametrische verdelingen komen, vanzelfsprekend, parameters voor en deze zullen moeten geschat worden. De meest effici¨ente methode om deze parameters te schatten is de maximumlikelihoodmethode. Om schattingen voor de parameters α, β, . . . te bekomen, maximaliseren we de (log-)likelihood l(α, β, . . . ; ~y ) = log Y fYi (yi ; α, β, . . .). Ten slotte moet het model gevalideerd worden. Het is dus belangrijk om de goodness-offit te testen, zo gaat men na hoe goed het model de gegeven observaties weergeeft. De schaarste aan gegevens, waarmee men vaak geconfronteerd wordt bij indexmodellering, kan het moeilijk maken om onderscheid te maken tussen een slechte fit door de steekproeffout enerzijds en een slechte fit die te wijten is aan de keuze van een verkeerd model anderzijds. Dit maakt het bijzonder belangrijk om in te schatten hoe veel variatie verwacht wordt door de steekproeffout en om diverse controles van het model toe te passen. We kunnen besluiten dat we de volgende stappen moeten doorlopen om de prijs van het temperatuurderivaat te bepalen: 66 Hoofdstuk 3. Temperatuurderivaten prijzen • De waarden van de index worden willekeurig getrokken uit een geschikte verdeling. • De payoffs worden bepaald aan de hand van de verkregen indices in de eerste stap. • Bepaal de gemiddelde payoff en verdisconteer dit resultaat. Tot slot geven we nog mee hoe Schiller et al. (2012) de aanpak via indexmodellering beschrijven. Zij stellen dat modellering van de index verder bouwt op de aanpak via historical burn analysis, door ook de verdeling van de temperatuurindex te schatten. Wanneer de verdeling relatief goed geschat kan worden, zal deze methode een stabieler resultaat opleveren. Wanneer we de drie onderstellingen uit de sectie 3.2 aanvullen met een vierde, namelijk dat de storingstermen εi , i = 1, . . . , n+1 onafhankelijk en identiek normaal verdeeld zijn, verkrijgen we een verdelingsfunctie voor de index Yi . Immers wanneer εi ∼ N (0, σ 2 ), dan geldt Yi ∼ N (β0 + β1 i, σ 2 ). Merk op dat de hypothese van een normale verdeling d n+1 − Ybn+1 ), d.i. de variantie het vaakst zal gemaakt worden. We kunnen nu ook Var(Y van de fout van de voorspelling (Yn+1 − Ybn+1 ), schatten9 . Zoals in Goetghebeur (2011) geldt ! n+1 2 n + 1 − 1 2 Var(Yn+1 − Ybn+1 ) = σ 2 · 1 + + Pn n+1 2 n i=1 i − 2 n+1 2 ( 2 ) 1 = σ 2 · 1 + + P n (n+1)2 n 2 i=1 i − (n + 1)i + 4 2 1 (n + 1) = σ 2 · 1 + + P (n+1)2 n 4 n 2 i=1 i − (n + 1)i + 4 2 1 (n + 1) 2 P P = σ · 1 + + Pn 2 n 4 i=1 i − 4(n + 1) ni=1 i + ni=1 (n + 1)2 ! 2 n+1 (n + 1) = σ2 · + (n+1)(2n+1)n n 4 − 4(n + 1) n(n+1) + n(n + 1)2 6 2 ! (n + 1) n+1 = σ2 · + (2n+1)n n 4 6 − 4 n(n+1) + n(n + 1) 2 ! n + 1 (n + 1) = σ2 · + n n 32 (2n + 1) − 2(n + 1) + n + 1 9 Er wordt opnieuw vertrokken van2 het simpel lineair regressiemodel Yi = β0 + β1 xi + εi , waaruit (X −¯ x) 2 b Var(Yp − Yp ) = σ · 1 + 1 + Pn p bekomen wordt. 2 n x) i=1 (xi −¯ 3.4. Dagelijkse simulatie 67 (n + 1) = σ2 · 4 n + 23 − (n + 1) 3 ! (n + 1) n + 1 + = σ2 · n n 13 (n − 1) n + 1 3(n + 1) 2 =σ · + n n(n − 1) (n − 1)(n + 1) + 3(n + 1) 2 =σ · n(n − 1) (n − 1 + 3)(n + 1) 2 =σ · n(n − 1) (n + 2)(n + 1) = σ2 · . n(n − 1) n+1 + n n ! We willen echter een schatting voor deze variantie, d.w.z. dat we de waarde σ 2 zullen n 1 X (Yi − Ybi )2 , dit moeten schatten. De kleinste kwadraten schatter van σ 2 is σ ˆ2 = n − 2 i=1 is tevens een onvertekende schatter. Door gebruik te maken van deze schatter bekomen we dus d n+1 − Ybn+1 ) = (n + 2)(n + 1) σ Var(Y ˆ2. n(n − 1) 3.4 Dagelijkse simulatie In deze sectie bespreken we een methode die nauwkeuriger is dan HBA en indexmodellering, we hebben ons hiervoor gebaseerd op Alexandridis and Zapranis (2013) en Schiller et al. (2012). Dagelijkse simulatie maakt gebruik van stochastische methoden om de temperatuur te modelleren op een dagelijkse basis. De temperatuur kan door zowel een discreet als een continu proces gemodelleerd worden. De dynamische modellen die hierdoor ontstaan, simuleren het toekomstig gedrag van de temperatuur rechtstreeks. De geschatte modellen kunnen uiteindelijk gebruikt worden om er de corresponderende indices uit af te leiden en om verschillende temperatuurderivaten te prijzen. Het effectief prijzen van de derivaten, steunend op dagelijkse simulatie, wordt besproken in Sectie 3.4.3 en 3.4.4. Dat men door het gebruik van modellen voor de dagtemperatuur een nauwkeuriger resultaat bekomt dan bij HBA (zie Sectie 3.2) voor de prijs van een derivaat, is vooral te verklaren doordat de dagelijkse modellering de beschikbare historische gegevens volledig gebruikt. Via dagelijkse modellering van de temperatuur kan eveneens een nauwkeurigere prijs bekomen worden dan via indexmodellering (zie Sectie 3.3). Zo gaat er 68 Hoofdstuk 3. Temperatuurderivaten prijzen bijvoorbeeld veel informatie (over zowel gewone als extreme gebeurtenissen) verloren wanneer de temperatuurindex berekend wordt via (bijvoorbeeld) een normaal verdeeld proces. Op de indices zijn er bijvoorbeeld bepaalde grenzen van toepassing, de degree days zijn begrensd door nul. Dagelijkse simulatie heeft nog een ander voordeel. Daar bij indexmodellering voor elke index steeds een ander model moet worden geschat, heeft men bij dagelijkse simulatie genoeg aan ´e´en model dat aan de data is aangepast. Dit model kan gebruikt worden voor alle beschikbare contracten op de markt op dezelfde locatie. Bovendien kunnen we via het gebruik van een dagelijks model een nauwkeurige representatie van alle indices en hun verdeling bekomen. Ten slotte is het, in tegenstelling tot indexmodellering en HBA, gemakkelijk om weersvoorspellingen te integreren. Het afleiden van een correct model voor de dagtemperatuur is echter geen eenvoudige opdracht (we verwijzen hiervoor ook naar Hoofdstuk 2). Observaties van de temperatuur vertonen seizoensgevoeligheid in het gemiddelde, de variantie, de verdeling en de autocorrelatie. Soms vertonen de waarnemingen ook een stijgende trend. Het is dan ook belangrijk dat het model deze eigenschappen correct weergeeft. Het risico dat aan de methode via dagelijkse simulatie verbonden is, is dat kleine fouten in het model kunnen leiden tot grote fouten in de prijzen van de contracten. In de literatuur worden twee methodes voorgesteld om de gemiddelde dagtemperatuur te modelleren, men maakt gebruik van ofwel een discreet ofwel een continu proces. 3.4.1 Discreet proces Men kan argumenteren dat het voordeliger is om een discreet proces te gebruiken om de temperatuur te modelleren. Men kan bijvoorbeeld stellen dat er onmiddellijk een discreet proces moet worden gebruikt, omdat de waarden van de temperatuur in discrete vorm gegeven worden. Een aanpak via discrete processen valt samen met een werkwijze via tijdreeksen. De algemene modellen die in dit kader gebruikt worden, werden besproken in Sectie 2.1. Alexandridis and Zapranis (2013) geven een overzicht van de modellen, gebaseerd op een discreet proces, die in de praktijk gehanteerd worden. Wij hebben er voor gekozen om het model voorgesteld door Sean D. Campbell en Francis X. Diebold uitgebreid te bespreken (zie Sectie 2.4.1), omdat dit model vaak geciteerd wordt in de literatuur. 3.4. Dagelijkse simulatie 3.4.2 69 Continu proces De continue modellen maken gebruik van stochastische differentiaalvergelijkingen, dit zijn vaak varianten van een mean-reverting Ornstein-Uhlenbeck stochastische differentiaalvergelijking. De temperatuur toont seizoensgevoeligheid in het gemiddelde en in de variantie, waardoor de component die de seizoensgevoeligheid beschrijft accuraat moet gemodelleerd worden. Dit werd ook reeds besproken in Sectie 2.3 en 2.4. Het is ook belangrijk dat een precieze schatting van de snelheid van mean-reversion verkregen wordt en dat de verdeling van de residuen correct geselecteerd wordt. Tot slot moet de geschikte lengte van de gebruikte tijdreeks van de temperatuur gekozen worden (zie hieromtrent ook de bespreking in Sectie 3.3) om de verschillende parameters van de gediscretiseerde versie van het model te kunnen schatten. Zodra het proces geschat is, kan men de waarde van een voorwaardelijke vordering bepalen door de verwachte waarde van de verdisconteerde toekomstige payoffs te beschouwen. Door de vaak complexe vorm van het proces en de pad-afhankelijke structuur van de (meeste) payoffs, heeft de uitdrukking voor de prijs van het contract doorgaans geen oplossingen die in gesloten vorm te schrijven zijn. In dat geval worden Monte Carlo-simulaties gebruikt. Dit omvat gewoonlijk het genereren van een groot aantal gesimuleerde scenario’s van de temperatuur(indices) om zo de mogelijke payoffs van het temperatuurderivaat te bepalen. De uiteindelijke prijs van het derivaat is dan het gemiddelde van alle gesimuleerde payoffs, eventueel verdisconteerd (om rekening te houden met de tijdwaarde van het geld). De precisie van de Monte Carlo-simulatie hangt natuurlijk af van de keuze van het temperatuurproces. In Alexandridis and Zapranis (2013) wordt een mooi overzicht gegeven van de in de literatuur reeds gebruikte continue processen. Wij hebben er hier twee van in detail besproken, namelijk het Alaton-model en het Benth-model (zie Sectie 2.4.2 en 2.4.3, respectievelijk). Om af te sluiten merken we nog op dat wij twee modellen gekozen hebben die gebruik maken van een Brownse beweging, er bestaan echter ook modellen die een gefractioneerde Brownse Beweging of een L´evy-proces hanteren. Het bouwen van een algoritme dat de basiskenmerken van de temperatuur correct zou defini¨eren, zou leiden tot een nauwkeurige prijsbepaling voor weerderivaten. In de volgende secties beschrijven we hoe we deze modellen nu ook concreet kunnen gebruiken om temperatuurderivaten te prijzen. 70 3.4.3 Hoofdstuk 3. Temperatuurderivaten prijzen Prijzen op basis van discrete dagelijkse simulatie In deze sectie zullen we aantonen hoe de discrete, econometrische modellen voor de temperatuur kunnen gebruikt worden om de prijs voor een temperatuurderivaat te bepalen. We zullen hiertoe het model voor de Belgische temperatuurgegevens, dat we in Sectie 2.4.4.2 opgesteld hebben, gebruiken. Aan de hand van dit model zullen we de prijs voor de calloptie, die we ook in Sectie 3.2.1 gebruikt hebben, bepalen. De calloptie is geschreven op de CAT-index over de periode van 1 mei tot en met 30 september. We veronderstellen dat het uitoefenniveau 2 500 is en dat de tick size e 20 is. De payoff van de optie wordt gegeven door (3.4). Hieronder zullen we beschrijven hoe we te werk gaan om de prijs van deze optie te bepalen, wanneer we het contract willen kopen voor de zomermaanden in 2011. We verwijzen naar Appendix B.4 voor de gebruikte R-code. Omdat we de prijs van de calloptie voor de zomermaanden in 2011 willen bepalen, hebben we een voorspelling voor de temperatuur in deze maanden nodig. Deze temperaturen kunnen we met behulp van het model simuleren. Hiertoe hebben we eerst 1 000 simulaties van het ARMA(2,2)-model, dat we gespecifieerd hadden voor de variabele restemp, voor de periode mei-september 2011 gemaakt. Om de correcte temperaturen te bekomen, hebben we vervolgens deze gesimuleerde waarden vermeerderd met de seizoenseffecten fs voor de periode mei-september. Ten slotte hebben we deze temperaturen bij elkaar opgeteld om de CAT-index te bekomen. We beschikken nu dus over 1 000 gesimuleerde waarden voor de CAT-index over de periode mei-september. Via deze CAT-index kunnen we de payoff van de calloptie (3.4) bepalen. We bezitten nu 1 000 mogelijke waarden voor de payoff van de calloptie en van deze waarden kunnen we uiteindelijk het gemiddelde nemen. Dit gemiddelde beschouwen we als de werkelijke payoff van het contract, we bekomen e 968.46. Wanneer we veronderstellen dat r = 2.25% de constante intrestvoet is en D(t, T ) = e−r(T −t) noteren, kunnen we de prijs voor de calloptie op de ingangsdatum van het contract bepalen. Voor de prijs van de calloptie op tijdstip t geldt P (t) = D(t, T ) · payoff, waarbij D(t, T ) de verdisconteringsfactor van maturiteit T naar tijdstip t is. Hierdoor is de prijs van de calloptie op de ingangsdatum van het contract gelijk aan 153 P (t) = e−0.0225 365 · 968.46 = 959.37. De prijs van de calloptie is e 959.37. Via HBA bekomen we voor hetzelfde contract e 845.98 (zie Sectie 3.2.1). De prijs voor de optie die we met behulp van een model voor 3.4. Dagelijkse simulatie 71 de Belgische temperatuur bekomen, is dus duurder dan de prijs die we verkrijgen via HBA. 3.4.4 Prijzen op basis van continue dagelijkse simulatie In deze sectie, gebaseerd op Alexandridis and Zapranis (2013) en Benth and Benth (2005), zullen we de formules voor het prijzen van verschillende temperatuurderivaten presenteren. We zullen hier de derivaten effectief gaan prijzen gebruik makend van dagelijkse simulatie aan de hand van een continu proces. Aangezien we ons hier toespitsen op de continue processen zullen we, wanneer we gebruik maken van de indices uit Hoofdstuk 1, de som in de definities vervangen door een integraal. Wanneer de markt compleet is, kan een unieke risiconeutrale kansmaat Q ∼ P verkregen worden, waarbij P de kansmaat in de echte wereld voorstelt. Deze verandering van maat maakt van het stochastisch proces een martingaal. Hierdoor kunnen financi¨ele derivaten geprijsd worden door de verdisconteerde verwachte waarde van de payoff van het derivaat onder de risiconeutrale maat te nemen. De markt voor het weer is echter incompleet, in de zin dat de onderliggende van het weerderivaat niet kan opgeslagen en verhandeld worden. Bovendien is de markt erg illiquide. In principe kan (uitgebreide) risiconeutrale waardering nog steeds uitgevoerd worden in incomplete markten. Echter, in incomplete markten kan een unieke prijs niet verkregen worden met behulp van de niet-arbitrage aanname. Met andere woorden, er bestaan veel equivalente martingalen en daardoor kunnen alleen boven- en ondergrenzen voor de prijzen van voorwaardelijke vorderingen worden ingevoerd. De verandering van maat van de echte wereld naar de risiconeutrale wereld onder de dynamiek van een Brownse beweging kan worden uitgevoerd door gebruik te maken van de stelling van Girsanov (cf. Shreve (2004)). De stelling van Girsanov transformeert een Brownse beweging onder P naar een nieuwe Brownse beweging onder Q en vertelt ons hoe een stochastisch proces verandert wanneer de maat wijzigt. Door het gebruik van een Brownse beweging veronderstellen we een normale verdeling, die behouden blijft door de verandering van maat. Het weer wordt echter niet steeds gedreven door een normale verdeling. Zo verkrijgt men vaak een beter model voor de temperatuur wanneer men uitgaat van een ingewikkelder sprongproces, zoals het L´evy-proces. Wanneer we gebruik maken van een L´evy-proces kan de verandering van maat worden uitgevoerd door gebruik te maken van de Esscher-transformatie, dit is een uitbreiding van de stelling van 72 Hoofdstuk 3. Temperatuurderivaten prijzen Girsanov10 . De Esscher-transformatie behoudt op haar beurt de eigenschappen van de verdeling van het L´evy-proces. De Esscher-transformatie wordt uitgebreider besproken in Sectie 3.4.4.2. Tot slot kan de geactualiseerde verwachte payoff van de verschillende contracten worden geschat. De oplossing van de stochastische differentiaalvergelijking die de dynamiek van de temperatuur beschrijft, moet echter bepaald worden om de verwachte payoff van elk derivaat te schatten. 3.4.4.1 Brownse beweging In deze sectie nemen we aan dat een normale verdeling gerechtvaardigd is voor de dynamiek van de temperatuur. 3.4.4.1.1 Stochastische differentiaalvergelijking voor de temperatuur Dagelijkse simulatie aan de hand van continue processen maakt gebruik van stochastische differentiaalvergelijkingen om de temperatuur te modelleren. Deze vergelijkingen zijn vaak varianten van een mean-reverting Ornstein-Uhlenbeck-proces. Het voorgestelde model kan worden gebruikt om de theoretische prijzen van de temperatuurderivaten te verkrijgen. Echter, om over te gaan tot de prijsstelling van de contracten, moet de oplossing van de stochastische differentiaalvergelijking bepaald worden. We volgen Alexandridis and Zapranis (2013) in hun keuze van een Ornstein-Uhlenbeckproces, met een tijdsafhankelijke snelheid van mean-reversion, om de temperatuur te modelleren. Door onze keuze van definities komt dit proces overeen met wat in deze masterproef het Hull-White-proces genoemd wordt. Stelling 2. Veronderstel dat de temperatuur een Hull-White-proces volgt dT (t) = dS(t) + κ(t)(T (t) − S(t))dt + σ(t)dW (t), met S(t), κ(t) en σ(t) deterministische functies. S(t) beschrijft de seizoenseffecten van de temperatuur, κ(t) staat voor de tijdsafhankelijke snelheid van mean-reversion en de volatiliteit σ(t) is een meetbare en begrensde functie. Een expliciete oplossing van deze vergelijking wordt gegeven door T (t) = S(t) + e Rt 0 κ(u)du Rt (T (0) − S(0)) + e 0 κ(u)du Z t σ(s)e− Rs 0 κ(u)du dW (s). 0 10 Men kan ook zeggen dat de stelling van Girsanov een speciaal geval is van de Esscher-transformatie voor de Brownse beweging. 3.4. Dagelijkse simulatie 73 Bewijs. Het bewijs van deze stelling werd reeds gegeven in Sectie 2.3.2, wanneer we X(t) = T (t) − S(t), α(t) = 0 en −β(t) = κ(t) stellen. Het resultaat van Stelling 2 zal rechtstreeks worden gebruikt bij de bepaling van de theoretische prijs van de temperatuurderivaten geschreven op CAT-, HDD-, CDD-, en PR-indices. 3.4.4.1.2 CAT- en PR-indices: futures en opties In deze sectie zullen we eerst een wiskundige uitdrukking opstellen voor de CAT-futureprijs. We weten reeds dat het niet mogelijk is om in incomplete markten een unieke prijs te verkrijgen met behulp van de niet-arbitrage aanname. Temperatuurderivaten worden geschreven op een temperatuurindex, dit is een actief dat niet kan worden opgeslagen of verhandeld. Om de formule voor het prijzen van derivaten te bepalen, zullen we eerst een risiconeutrale kansmaat Q ∼ P moeten vinden, waaronder alle activa na verdiscontering martingalen zijn. In het geval van weerderivaten is elke equivalente11 maat Q een risiconeutrale kansmaat. Als Q de risiconeutrale kansmaat en r de constante samengestelde intrestvoet is, dan wordt de niet-arbitrage futureprijs van een CAT-contract met waarde Z τ2 T (τ )dτ op tijdstip τ2 gegeven door (cf. Shreve (2004)) τ1 Z τ2 FCAT (t, τ1 , τ2 ) = EQ τ1 T (τ )dτ Ft , t ≤ τ1 < τ2 . f Via Z t de stelling van Girsanov hebben we, onder de equivalente maat Q, W (t) = W (t) − θ(u)du of nog 0 f (t) = dW (t) − θ(t)dt, dW (3.7) waar θ(t) een meetbare en begrensde re¨eelwaardige functie is, die de marktprijs voor risico voorstelt. Aangezien temperatuur niet verhandelbaar is, moet de marktprijs voor risico worden opgenomen in het waarderingsmodel. De keuze van θ(t) geeft eigenlijk weer wat de voorkeur voor risico van de spelers in de markt is. Het weerspiegelt, met andere woorden, de mening van de handelaar over het zich blootstellen aan risico’s. Zoals blijkt uit (3.7) is de verandering van maat van een stochastisch proces van een actief nauw verwant aan het concept marktprijs voor risico. Uit (3.8) hieronder kunnen we 11 Zij Ω een niet-ledige verzameling en F een σ-algebra van deelverzamelingen van Ω. Twee kansmaten P en Q op (Ω,F) zijn equivalent als ze dezelfde verzamelingen in F een kans nul toemeten (cf. Shreve (2004)). 74 Hoofdstuk 3. Temperatuurderivaten prijzen afleiden dat de drift-snelheid van het stochastisch proces van de temperatuur gecorrigeerd wordt met een parameter die de marktprijs voor risico reflecteert. Tot nog toe werd in de meeste studies aangenomen dat de marktprijs voor risico gelijk is aan nul. Onlangs deden velen echter onderzoek naar de marktprijs voor risico en men vond dat θ(t) niet steeds gelijk is aan nul. In Alexandridis and Zapranis (2013) wordt een mooi overzicht gegeven van de studies en de methoden die gebruikt worden om de marktprijs voor risico te schatten. De meest voorkomende aanpak is deze die Alaton et al. (2002) voorstellen. Zij suggereren dat de marktprijs voor risico kan worden geschat op basis van de historische marktgegevens. Meer specifiek kan θ(t) worden berekend door te kijken naar de marktprijs van de contracten. De marktprijs voor risico is de waarde die maakt dat de prijs van het theoretisch model gelijk is aan de waarneembare marktprijs. Als we (3.7) substitueren in het Hull-White-proces voor de temperatuur, krijgen we het stochastisch proces voor de temperatuur onder de risiconeutrale kansmaat Qθ f (t) + θ(t)dt) dT (t) = dS(t) + κ(t)(T (t) − S(t))dt + σ(t)(dW (3.8) f (t). = dS(t) + [κ(t)(T (t) − S(t)) + σ(t)θ(t)]dt + σ(t)dW Merk op dat Q de risiconeutrale kansmaat is, waarbij Q ∼ P. Omdat we in deze sectie gebruik maken van een Brownse beweging staat Qθ voor een subklasse van deze kansen bepaald door de stelling van Girsanov. Aangezien we enkel in deze kansen ge¨ınteresseerd zijn, zullen we om de notatie te vereenvoudigen deze subklasse van kansmaten ook met Q noteren. De oplossing van de voorgaande stochastische differentiaalvergelijking onder Q is Z t Rt Rt Rs κ(u)du κ(u)du 0 0 T (t) =S(t) + e (T (0) − S(0)) + e σ(s)θ(s)e− 0 κ(u)du ds 0 (3.9) Z t Rs Rt − 0 κ(u)du f κ(u)du 0 +e σ(s)e dW (s). 0 Het bewijs van deze bewering werd reeds gegeven in Sectie 2.3.2, wanneer we X(t) = T (t) − S(t), α(t) = σ(t)θ(t) en −β(t) = κ(t) stellen. Gebruik makend van deze laatste uitdrukking, vinden we de prijs van een futurecontract op de CAT-index op tijdstip t, waarbij t ≤ τ1 < τ2 . Stelling 3. De CAT-futureprijs voor t ≤ τ1 < τ2 wordt gegeven door Z τ2 Z τ 2 FCAT (t, τ1 , τ2 ) = EQ T (s)ds Ft = S(s)ds + I1 + I2 , τ1 waarbij Z τ2 I1 = Te(t) e τ1 Rs t κ(z)dz ds, τ1 (3.10) 3.4. Dagelijkse simulatie τ1 Z Z τ2 Rs e I2 = t 0 κ(z)dz 75 σ(u)θ(u)e R0 u κ(z)dz Z τ2 τ2 Z 0 κ(z)dz σ(u)θ(u)e R0 u κ(z)dz dsdu, u τ1 τ1 Rs e dsdu + met Te(t) = T (t) − S(t). Bewijs. Per definitie hebben we Z τ2 FCAT (t, τ1 , τ2 ) = EQ Zτ1τ2 T (s)ds Ft (S(s) + T (s) − S(s))ds Ft = EQ τ Z τ 2 Z τ2 1 e S(s)ds + EQ T (s)ds Ft , = τ1 τ1 waarbij Te(s) = T (s)−S(s). Uit de oplossing van de stochastische differentiaalvergelijking onder Q (3.9), kunnen we afleiden dat er voor t ≤ s geldt Te(s) = T (s) − S(s) Rs Rs Z s Ru (T (t) − S(t)) + e t σ(u)θ(u)e− t κ(z)dz du t Z s Ru Rs f (u) σ(u)e− t κ(z)dz dW + e t κ(z)dz t Z s Z s Rs Rs Rs κ(z)dz κ(z)dz e f (u). u t σ(u)e u κ(z)dz dW σ(u)θ(u)e du + =e T (t) + =e t κ(z)dz κ(z)dz t t Verder geldt er dat i h EQ Te(s) Ft R Z s Z s Rs Rs s κ(z)dz κ(z)dz κ(z)dz f (u) Ft = EQ e t Te(t) + σ(u)θ(u)e u du + σ(u)e u dW t t Z s Rs Rs = e t κ(z)dz Te(t) + σ(u)θ(u)e u κ(z)dz du. t Voor de CAT-futureprijs geldt FCAT (t, τ1 , τ2 ) Z τ 2 Z τ2 e = S(s)ds + EQ T (s)ds Ft Zτ1τ2 Z τ2 τ1 h i = S(s)ds + EQ Te(s) Ft ds τ τ Z 1τ2 Z 1τ2 R Z s Rs s κ(z)dz κ(z)dz = S(s)ds + et Te(t) + σ(u)θ(u)e u du ds τ1 τ1 t Z τ2 Z τ2 R Z τ2 Z s Rs s κ(z)dz e t = S(s)ds + e T (t)ds + σ(u)θ(u)e u κ(z)dz duds τ1 τ1 τ1 t 76 Hoofdstuk 3. Temperatuurderivaten prijzen Z τ2 S(s)ds + I1 + I2 . = τ1 Z τ2 Rs e Er geldt dus I1 = Te(t) t κ(z)dz Z τ2 s Z Rs σ(u)θ(u)e ds en I2 = κ(z)dz duds. Enkel t τ1 τ1 u de uitdrukking voor I2 moet nog herschreven worden. Wanneer we de karakteristieke functie I defini¨eren ID (x) = ( 1 x∈D 0 x∈ /D kunnen we I2 als volgt omvormen Z τ2 Z s Rs σ(u)θ(u)e u κ(z)dz duds I2 = Zτ1τ2 Zt τ2 Rs I[t,s] (u)σ(u)θ(u)e u κ(z)dz duds = Zτ1τ2 Zt τ2 Rs = I[t,s] (u)σ(u)θ(u)e u κ(z)dz dsdu Zt τ1 Zτ1τ2 Z Rs κ(z)dz = I[t,s] (u)σ(u)θ(u)e u dsdu + t τ1 τ2 τ1 Z , τ2 I[t,s] (u)σ(u)θ(u)e Rs u κ(z)dz dsdu. τ1 De tweede term in deze uitdrukking is nul als s kleiner is dan u, daardoor kunnen we de grenzen van de binnenste integraal aanpassen Z τ1 Z τ2 Z τ 2 Z τ2 Rs Rs κ(z)dz u I2 = σ(u)θ(u)e u κ(z)dz dsdu σ(u)θ(u)e dsdu + τ1 u Zt τ1 Zτ1τ2 R Z τ2 Z τ2 R R0 R s s 0 κ(z)dz κ(z)dz 0 u = e 0 κ(z)dz σ(u)θ(u)e u κ(z)dz dsdu, e σ(u)θ(u)e dsdu + t τ1 τ1 u en dit is de gezochte uitdrukking voor I2 . Deze stelling geeft de CAT-futureprijs op tijdstip t ≤ τ1 < τ2 . Met andere woorden, de prijs van een CAT-future v´oo´r de contractperiode. Hierdoor is (3.10) een waardering buiten de contractperiode. Teneinde de futureprijs te evalueren binnen de contractperiode kan de bovenstaande formule eenvoudig aangepast worden. Stelling 4. De CAT-futureprijs voor τ1 ≤ t ≤ τ2 wordt gegeven door Z t FCAT (t, τ1 , τ2 ) = T (s)ds + FCAT (t, t, τ2 ). τ1 Bewijs. Voor de CAT-futureprijs geldt Z FCAT (t, τ1 , τ2 ) = EQ τ2 τ1 T (s)ds Ft 3.4. Dagelijkse simulatie 77 t τ2 T (s)ds Ft = EQ T (s)ds + t τ1 Z τ 2 Z t T (s)ds Ft = T (s)ds + EQ t τ Z 1t T (s)ds + FCAT (t, t, τ2 ), = Z Z τ1 Z t T (s)ds Ft -meetbaar is. de voorlaatste gelijkheid geldt omdat τ1 In de volgende stelling bekijken we de dynamiek van de CAT-futureprijs onder de risiconeutrale kansmaat Q. Stelling 5. De dynamiek van FCAT (t, τ1 , τ2 ) onder de kansmaat Q is f (t), dFCAT (t, τ1 , τ2 ) = ΣCAT (t, τ1 , τ2 )dW Z τ2 R s waarbij ΣCAT (t, τ1 , τ2 ) = σ(t) e t κ(z)dz ds. τ1 Bewijs. De futureprijs FCAT (t, τ1 , τ2 ) is een martingaal onder de risiconeutrale maat Q (cf. Shreve (2004)), hierdoor weten we dat er in de uitdrukking voor dFCAT (t, τ1 , τ2 ) geen dt-term zal staan. Uit Stelling 3 weten we dat het volgende geldt Z τ2 Z τ2 dFCAT (t, τ1 , τ2 ) = d S(s)ds + I1 + I2 = d S(s)ds + dI1 + dI2 . τ1 τ1 De eerste term in het rechterlid is gelijk aan nul. We maken gebruik van de formule van Itˆo om dI1 te bepalen dI1 Z e = d T (t) τ2 Rs e t κ(z)dz ds τ1 = d(Te(t)f (t)) = Te(t)df (t) + f (t)dTe(t) + dTe(t)df (t) Z τ2 R Z τ2 R Z s s κ(z)dz κ(z)dz e e e t t = T (t) e (−κ(t)dt)ds + e dsdT (t) + dT (t) τ1 τ1 τ2 e Rs t κ(z)dz (−κ(t)dt)ds. τ1 Omdat er in de uitdrukking voor dFCAT (t, τ1 , τ2 ) geen dt-term zal staan, kunnen we dit ook als volgt noteren Z τ2 dI1 = . . . dt + Rs e t κ(z)dz dsdTe(t). τ1 De uitdrukking voor dI2 zal enkel bestaan uit een dt-term, we noteren dI2 = . . . dt. 78 Hoofdstuk 3. Temperatuurderivaten prijzen Hierdoor krijgen we τ2 Z e dFCAT (t, τ1 , τ2 ) = . . . dt + Rs t κ(z)dz dsdTe(t), τ1 uit (3.8) volgt dan Z τ2 e dFCAT (t, τ1 , τ2 ) = . . . dt + Zτ1τ2 e = . . . dt + Rs t Rs t κ(z)dz f (t) ds . . . dt + σ(t)dW κ(z)dz f (t). dsσ(t)dW τ1 Omdat FCAT (t, τ1 , τ2 ) een martingaal is onder de risiconeutrale maat volgt het gestelde. Nu kunnen we van de voorgaande stelling gebruik maken om de prijs van een calloptie geschreven op CAT-futures te bepalen. Stelling 6. De prijs op tijdstip t ≤ τ van een calloptie geschreven op een CAT-future met uitoefenprijs K op uitvoeringstijdstip τ ≤ τ1 is CCAT (t, τ, τ1 , τ2 ) = e−r(τ −t) [(FCAT (t, τ1 , τ2 ) − K)Φ(d(t, τ, τ1 , τ2 )) + Φ0 (d(t, τ, τ1 , τ2 ))Σt,τ ] , waarbij FCAT (t, τ1 , τ2 ) − K , Σt,τ Z τ Σ2CAT (s, τ1 , τ2 )ds, = d(t, τ, τ1 , τ2 ) = Σ2t,τ t en Φ de cumulatieve standaard normale verdelingsfunctie is. Bewijs. De optieprijs wordt bij definitie gegeven via de risiconeutrale prijsformule CCAT (t, τ, τ1 , τ2 ) = e−r(τ −t) EQ [max{FCAT (τ, τ1 , τ2 ) − K, 0}| Ft ] . f (t) en via inteUit de vorige stelling weten we dat dFCAT (t, τ1 , τ2 ) = ΣCAT (t, τ1 , τ2 )dW gratie over het interval [t, τ ] kunnen we de dynamiek van de futureprijs als volgt neerschrijven Z τ FCAT (τ, τ1 , τ2 ) = FCAT (t, τ1 , τ2 ) + f (s). ΣCAT (s, τ1 , τ2 )dW t Hieruit blijkt dat FCAT (τ, τ1 , τ2 ) geconditioneerd op Z FCAT (t, τ1 , τ2 ) een normale verdeling τ Σ2CAT (s, τ1 , τ2 )ds. volgt met gemiddelde FCAT (t, τ1 , τ2 ) en variantie t Om de optieprijs expliciet te bepalen, gaan we te werk zoals in Kaas et al. (2009). Stel 3.4. Dagelijkse simulatie 79 dat X normaal verdeeld is met gemiddelde µ en variantie σ 2 , X ∼ N (µ, σ 2 ). Als U standaard normaal verdeeld is, U ∼ N (0, 1), dan geldt X = σU + µ ∼ N (µ, σ 2 ). We kunnen nu E[max{X − d, 0}] als volgt herschrijven E[max{X − d, 0}] = E[max{σU + µ − d, 0}] d−µ = σE max U − ,0 . σ Aangezien U standaard normaal verdeeld is, kunnen we het volgende neerschrijven Z ∞ max{u − t, 0}φ(u)du E[max{U − t, 0}] = −∞ Z ∞ (u − t)φ(u)du = Zt ∞ Z ∞ uφ(u)du − t φ(u)du, = t t waarbij φ(u) staat normale kansdichtheidsfunctie. We weten ook dat Z ∞ voor de standaard Z ∞ 0 de gelijkheid uφ(u)du = (−φ (u))du = φ(t) geldt, zodat er onmiddellijk volgt t t Z ∞ Z uφ(u)du − t E[max{U − t, 0}] = ∞ φ(u)du t t = φ(t) − t(1 − Φ(t)). Uit het bovenstaande kunnen we dan halen d−µ E[max{X − d, 0}] = σE max U − ,0 σ d−µ d−µ d−µ =σ φ − 1−Φ σ σ σ d−µ d−µ = σφ − (d − µ) 1 − Φ . σ σ We passen nu deze formule toe met de gegevens Z τ uit het bewijs, met andere woorden 2 X = FCAT (τ, τ1 , τ2 ), µ = FCAT (t, τ1 , τ2 ), σ = Σ2CAT (s, τ1 , τ2 )ds en d = K. Omdat X ∼ N (µ, σ 2 ) conditioneel op Ft , krijgen we t EQ [max{FCAT (τ, τ1 , τ2 ) − K, 0}| Ft ] K −µ K −µ = σφ − (K − µ) 1 − Φ σ σ K − FCAT (t, τ1 , τ2 ) K − FCAT (t, τ1 , τ2 ) = Σt,τ φ − (K − FCAT (t, τ1 , τ2 )) 1 − Φ . Σt,τ Σt,τ 80 Hoofdstuk 3. Temperatuurderivaten prijzen Wanneer we deze gelijkheid substitueren in de uitdrukking voor CCAT (t, τ, τ1 , τ2 ), rekening houdend met Φ0 (x) = φ(x) = φ(−x) en 1 − Φ(x) = Φ(−x), dan bekomen we de gezochte uitdrukking voor de optieprijs. We bekijken tot slot nog de future- en optieprijs voor een contract geschreven op de Pacific Rim index. Uit de definities volgt dat de Pacific Rim index gelijk is aan het gemiddelde van de CAT-index over een specifieke tijdsperiode. De prijs voor de PRfuture op tijdstip t ≤ τ1 < τ2 is per definitie (cf. Shreve (2004)) Z τ2 1 FPR (t, τ1 , τ2 ) = EQ T (s)ds Ft . τ 2 − τ 1 τ1 Hieruit kunnen we afleiden dat de volgende relatie tussen de CAT-futureprijs en de PR-futureprijs geldt 1 FCAT (t, τ1 , τ2 ). τ2 − τ1 We kunnen nu ook de prijs van een calloptie geschreven op een PR-future bepalen. Via FPR (t, τ1 , τ2 ) = de risiconeutrale prijsformule bekomen we de volgende relatie CPR (t, τ, τ1 , τ2 ) = e−r(τ −t) EQ [max{FPR (τ, τ1 , τ2 ) − K, 0}| Ft ] 1 −r(τ −t) =e EQ max FCAT (τ, τ1 , τ2 ) − K, 0 Ft τ2 − τ1 1 = e−r(τ −t) EQ [max{FCAT (τ, τ1 , τ2 ) − K(τ2 − τ1 ), 0}| Ft ] . τ2 − τ1 We kunnen nu het bewijs van Stelling 6 herhalen en komen zo tot de volgende uitdrukking voor de prijs van de calloptie CPR (t, τ, τ1 , τ2 ) 1 = e−r(τ −t) [(FCAT (t, τ1 , τ2 ) − K(τ2 − τ1 ))Φ(d(t, τ, τ1 , τ2 )) + Φ0 (d(t, τ, τ1 , τ2 ))Σt,τ ] , τ2 − τ1 met d(t, τ, τ1 , τ2 ) = 3.4.4.1.3 FCAT (t, τ1 , τ2 ) − K(τ2 − τ1 ) . Σt,τ HDD- en CDD-indices: futures en opties In deze sectie bespreken we de prijsformules voor futures en opties geschreven op HDDen CDD-indices. We presenteren analoge stellingen als in de vorige sectie. Als er een contract geschreven wordt op HDD- of CDD-indices, dan is dit voor een welbepaalde periode waarover de indices geaccumuleerd worden. De AccHDD- en de AccCDD-indices 3.4. Dagelijkse simulatie 81 over een periode [τ1 , τ2 ] worden gegeven door (in deze sectie noteren we voor de eenvoud HDD en CDD waar we AccHDD en AccCDD bedoelen) Z τ2 max{c − T (s), 0}ds, HDD = τ1 Z τ2 max{T (s) − c, 0}ds. CDD = τ1 Wij veronderstellen voor de referentietemperatuur c = 18◦ C, in de stellingen zullen we om deze zo algemeen mogelijk te houden echter c blijven noteren. Uit de definitie van HDD en CDD blijkt dat de prijsformules voor beide indices zeer gelijkaardig zullen zijn. We beginnen met het geven van een wiskundige uitdrukking voor de HDD-futureprijs. Als Q de risiconeutrale kansmaat en r de constante samengestelde intrestvoet Z τ2 is, dan max{c − wordt de niet-arbitrage futureprijs van een HDD-contract met waarde τ1 T (τ ), 0}dτ op tijdstip τ2 gegeven door (cf. Shreve (2004)) Z τ 2 FHDD (t, τ1 , τ2 ) = EQ max{c − T (τ ), 0}dτ Ft , t ≤ τ1 < τ2 . (3.11) τ1 Analoog wordt de niet-arbitrage futureprijs van een CDD-contract op tijdstip t ≤ τ1 < τ2 gegeven door Z τ2 FCDD (t, τ1 , τ2 ) = EQ τ1 max{T (τ ) − c, 0}dτ Ft . In de volgende stelling bespreken we de relatie tussen FCAT , FHDD en FCDD . Stelling 7. De CDD-, HDD- en CAT-futureprijzen zijn op de volgende manier aan elkaar gelinkt FHDD (t, τ1 , τ2 ) = c(τ2 − τ1 ) − FCAT (t, τ1 , τ2 ) + FCDD (t, τ1 , τ2 ). Bewijs. De HDD-futureprijs wordt gegeven door (3.11). In deze uitdrukking kunnen we max{c − T (τ ), 0} als volgt herschrijven max{c − T (τ ), 0} = c − T (τ ) + max{T (τ ) − c, 0}. En dus FHDD (t, τ1 , τ2 ) Z τ2 = EQ (c − T (τ ) + max{T (τ ) − c, 0}) dτ Ft Zτ1τ2 Z τ 2 Z = EQ cdτ Ft − EQ T (τ )dτ Ft + EQ τ1 τ1 = c(τ2 − τ1 ) − FCAT (t, τ1 , τ2 ) + FCDD (t, τ1 , τ2 ), dit is de gezochte uitdrukking. τ2 τ1 max{T (τ ) − c, 0}dτ Ft 82 Hoofdstuk 3. Temperatuurderivaten prijzen Deze stelling geeft aan dat de prijsformules van futures op CDD- en HDD-indices gelijkaardig en aan elkaar gelinkt zijn. We zullen hierdoor in de volgende stellingen enkel focussen op de prijsformules van contracten geschreven op de CDD-index. Stelling 8. De CDD-futureprijs voor 0 ≤ t ≤ τ1 < τ2 wordt gegeven door Z τ2 max{T (s) − c, 0}ds Ft FCDD (t, τ1 , τ2 ) = EQ τ1 Rs κ(z)dz e Z τ2 t m t, s, e T (t) ds, v(t, s)Ψ = v(t, s) τ1 (3.12) waarbij Z s Ru Rs Rs Rs κ(z)dz κ(z)dz κ(z)dz σ(u)θ(u)e− t κ(z)dz du − c, m t, s, e t Te(t) = S(s) + e t Te(t) + e t t Z s Rs Ru v 2 (t, s) = e2 t κ(z)dz σ 2 (u)e−2 t κ(z)dz du, t Te(t) = T (t) − S(t) en Ψ(x) = xΦ(x) + Φ0 (x) met Φ de cumulatieve standaard normale verdelingsfunctie. Bewijs. Uit (3.9) weten we dat T (s) onder de risiconeutrale kansmaat Q normaal verdeeld is met gemiddelde en variantie respectievelijk Rs EQ [T (s)| Ft ] = S(s) + e t κ(z)dz Te(t) + e Rs t κ(z)dz Z s σ(u)θ(u)e− Ru t κ(z)dz du, t 2 VarQ [T (s)| Ft ] = e Rs t κ(z)dz s Z σ 2 (u)e−2 Ru t κ(z)dz du. t Rs Hieruit volgt dat T (s) − c normaal verdeeld is met gemiddelde m t, s, e t κ(z)dz Te(t) en variantie v 2 (t, s). Nu kunnen we de CDD-futureprijs bepalen, er geldt Z τ 2 FCDD (t, τ1 , τ2 ) = EQ max{T (s) − c, 0}ds Ft Z τ2 τ1 = EQ [max{T (s) − c, 0}| Ft ] ds. τ1 We kunnen EQ [max{T (s) − c, 0}|Ft ] zoals in het bewijs van Stelling 6 bepalen EQ [max{T (s) − c, 0}|Ft ] c−µ c−µ = σφ − (c − µ) 1 − Φ σ σ Rs −m t, s, e t κ(z)dz Te(t) = v(t, s)φ v(t, s) 3.4. Dagelijkse simulatie 83 Rs −m t, s, e t κ(z)dz Te(t) + m t, s, e t κ(z)dz Te(t) 1 − Φ v(t, s) Rs Rs κ(z)dz e t m t, s, e t κ(z)dz Te(t) m t, s, e T (t) Rs + m t, s, e t κ(z)dz Te(t) Φ = v(t, s)Φ0 v(t, s) v(t, s) Rs m t, s, e t κ(z)dz Te(t) = v(t, s)Φ0 v(t, s) Rs Rs κ(z)dz e κ(z)dz e t t m t, s, e T (t) m t, s, e T (t) + v(t, s) Φ v(t, s) v(t, s) Rs κ(z)dz e t m t, s, e T (t) . = v(t, s)Ψ v(t, s) Rs Hieruit volgt het gestelde. De voorgaande stelling geeft de CDD-futureprijs op tijdstip t ≤ τ1 < τ2 , met andere woorden de prijs van een CDD-future v´o´or de contractperiode. Hierdoor bekomen we een waardering buiten de contractperiode. Teneinde de futureprijs te evalueren binnen de contractperiode kan de bovenstaande formule eenvoudig aangepast worden. Het bewijs van onderstaande stelling is analoog aan dat van Stelling 4. Stelling 9. De CDD-futureprijs voor τ1 ≤ t < τ2 wordt gegeven door Z t FCDD (t, τ1 , τ2 ) = max{T (s) − c, 0}ds + FCDD (t, t, τ2 ). τ1 Zoals voor CAT-contracten is het ook voor de CDD-futureprijs mogelijk om de dynamiek onder de risiconeutrale maat Q te bepalen, dit wordt in Stelling 10 besproken. Stelling 10. De dynamiek onder Q van FCDD (t, τ1 , τ2 ) voor 0 ≤ t ≤ τ1 < τ2 wordt gegeven door f (t), dFCDD (t, τ1 , τ2 ) = ΣCDD (t, τ1 , τ2 )dW waarbij Z τ2 ΣCDD (t, τ1 , τ2 ) = σ(t) τ1 Rs m t, s, e t κ(z)dz Te(t) ds, e t κ(z)dz Φ v(t, s) Rs met Φ de cumulatieve standaard normale verdelingsfunctie. 84 Hoofdstuk 3. Temperatuurderivaten prijzen Bewijs. De futureprijs FCDD (t, τ1 , τ2 ) is een Q-martingaal (cf. Shreve (2004)), hierdoor weten we dat er in de uitdrukking voor dF (t, τ1 , τ2 ) geen dt-term zal staan. Als we Rs !CDD Z τ2 κ(z)dz x m t, s, e t v(t, s)Ψ ds stellen, krijgen we via de formule van Itˆo f (t, x) = v(t, s) τ1 ∂f ∂f 1 ∂ 2f (t, X(t))dt + (t, X(t))dX(t) + (t, X(t))dX(t)dX(t). ∂t ∂x 2 ∂x2 Wanneer we X(t) = Te(t) stellen, volgt uit Stelling 8 dat we bovenstaande gelijkheid als df (t, X(t)) = volgt kunnen herschrijven dFCDD (t, τ1 , τ2 ) Z τ2 = . . . dt + τ1 Rs Rs m t, s, e t κ(z)dz Te(t) m0 t, s, e t κ(z)dz Te(t) v(t, s)Ψ0 dsdTe(t) v(t, s) v(t, s) 1 ∂ 2f (t, Te(t))dTe(t)dTe(t) + 2 ∂x2 Rs Z τ2 m t, s, e t κ(z)dz Te(t) Rs e t κ(z)dz dsdTe(t), = . . . dt + Φ v(t, s) τ1 f (t). In de laatste gelijkheid hebben we gebruik gemaakt met dTe(t) = . . . dt + σ(t)dW Rs Rs 0 0 κ(z)dz e 2 e e t van Ψ (x) = Φ(x), dT (t)dT (t) = σ (t)dt en m t, s, e T (t) = e t κ(z)dz , zowel Rs Ψ(x) als m t, s, e t κ(z)dz Te(t) werden gedefinieerd in Stelling 8. We tonen nog aan dat Ψ0 (x) = Φ(x) geldt. Per definitie is Ψ(x) gelijk aan xΦ(x) + Φ0 (x), er geldt dan Ψ0 (x) = Φ(x) + xΦ0 (x) + Φ00 (x) = Φ(x) + xφ(x) + φ0 (x) = Φ(x) + xφ(x) − xφ(x) = Φ(x), waarbij φ(x) de standaard normale kansdichtheidsfunctie is. f -term verschillend van nul is, geldt dus Omdat enkel de dW Rs Z τ2 R e(t) t κ(z)dz T m t, s, e s f (t). dsdW dFCDD (t, τ1 , τ2 ) = σ(t) e t κ(z)dz Φ v(t, s) τ1 Dit is de uitdrukking die we zochten. In de voorgaande stelling geeft de term ΣCDD (t, τ1 , τ2 ) de termijnstructuur van de volatiliteit van CDD-futures weer. Door hier gebruik van te maken, kunnen we de prijs van een calloptie op een CDD-future bepalen. 3.4. Dagelijkse simulatie 85 Stelling 11. De prijs op tijdstip t ≤ τ van een calloptie geschreven op een CDD-future met uitoefenprijs K en uitoefentijdstip τ ≤ τ1 wordt gegeven door Z τ2 −r(τ −t) CCDD (t, τ, τ1 , τ2 ) = e EQ max v(τ, s)Z t, s, τ, Te(t) ds − K, 0 , τ1 waarbij Z Rs κ(z)dz e e e t Z t, s, τ, T (t) = Ψ τ, s, e T (t) + τ Rs σ(u)θ(u)e u κ(z)dz du + Σ(s, t, τ )Y t met m(t, s, x) e Ψ(t, s, x) = Ψ , v(t, s) Z τ Rs 2 Σ (s, t, τ ) = σ 2 (u)e2 u κ(z)dz du, t en Y een standaard normaal verdeelde stochastische variabele. Bewijs. De optieprijs wordt via de risiconeutrale prijsformule gegeven door CCDD (t, τ, τ1 , τ2 ) = e−r(τ −t) EQ [max {FCDD (τ, τ1 , τ2 ) − K, 0}| Ft ] . e geldt Gelet op (3.12) en de notatie Ψ Z τ2 R e τ, s, e τs κ(z)dz Te(τ ) ds. FCDD (τ, τ1 , τ2 ) = v(τ, s)Ψ τ1 Uit (3.9) volgt dat we Te(τ ) als volgt kunnen noteren Te(τ ) Z τ Z τ Rτ R Rτ Ru κ(z)dz − tu κ(z)dz κ(z)dz e f (u) t t =e T (t) + e σ(u)θ(u)e du + e σ(u)e− t κ(z)dz dW t t Z τ Z τ Rτ Rτ Rτ f (u). σ(u)θ(u)e u κ(z)dz du + σ(u)e u κ(z)dz dW = e t κ(z)dz Te(t) + Rτ t κ(z)dz t t Wanneer we dit substitueren in de uitdrukking voor FCDD , krijgen we FCDD (τ, τ1 , τ2 ) Z τ2 Rs κ(z)dz e e τ = v(τ, s)Ψ τ, s, e T (τ ) ds τ1 Z τ2 Z Rs κ(z)dz e e = v(τ, s)Ψ τ, s, e t T (t) + τ1 t τ Rs σ(u)θ(u)e u κ(z)dz Z du + τ σ(u)e Rs u κ(z)dz t Dit invullen in de risiconeutrale prijsformule levert het gewenste resultaat. f dW (u) ds. 86 3.4.4.2 Hoofdstuk 3. Temperatuurderivaten prijzen L´ evy-proces In deze sectie vertrekken we van het idee dat de temperatuur gemodelleerd wordt door een Hull-White-proces dat gedreven wordt door een L´evy-proces in plaats van een Brownse beweging. We geven hiertoe eerst de definitie van een L´evy-proces zoals in Schoutens (2003). Veronderstel dat φ de karakteristieke functie van een verdeling is, deze functie bepaalt de verdelingsfunctie F op een unieke manier. Definitie 13. De karakteristieke functie φ van een verdeling, of equivalent een stochastische variabele X, is de Fourier-Stieltjes transformatie van de verdelingsfunctie F (x) = P(X ≤ x) Z +∞ φX (u) = E[exp(iuX)] = exp(iux)dF (x). −∞ Als, voor elk positief geheel getal n, φ(u) ook de n-de macht van een karakteristieke functie is, dan zeggen we dat de verdeling oneindig deelbaar is. Voor zulke oneindig deelbare verdelingen kunnen we een stochastisch proces X = {Xt , t ≥ 0} defini¨eren, met name het L´evy-proces. Definitie 14. Een L´ evy-proces is een stochastisch proces X dat start in nul en onafhankelijke en stationaire toenames heeft zodanig dat de verdeling van een toename over [s, s + t] met s, t ≥ 0, i.e. Xt+s − Xs , (φ(u))t als karakteristieke functie heeft. Definitie 15. Beschouw een re¨eelwaardige functie f : [a, b] → R. Veronderstel dat voor elke t ∈ (a, b] de functie f rechtscontinu is en een linkerlimiet heeft, dan zeggen we dat de functie c` adl` ag is, van het Franse continue `a droite et limites `a gauche. Elk L´evy-proces heeft een c`adl`ag-modificatie die op haar beurt een L´evy-proces is. Het is standaard om met deze c`adl`ag-versie van het proces te werken. Op die manier zijn de gesimuleerde paden van een L´evy-proces bijna zeker continu van rechts en hebben ze limieten van links. We merken ook nog op dat de cumulant karakteristieke functie ψ(u) = ln(φ(u)) vaak de karakteristieke exponent wordt genoemd en voldoet aan de L´evy-Khintchine-formule Z ∞ 1 2 2 ψ(u) = iγu − σ u + (exp(iux) − 1 − iuxI|x|<1 )ν(dx), 2 −∞ met γ ∈ R, σ 2 ≥ 0 en ν de L´evy-maat. 3.4. Dagelijkse simulatie 87 Definitie 16. Zij X = {Xt , t ≥ 0} een L´evy-proces en zij ν een maat op R\{0} met Z ∞ inf{1, x2 }ν(dx) < ∞, −∞ dan wordt ν de L´ evy-maat genoemd. We zeggen dat de oneindig deelbare verdeling een L´evy-triplet [γ, σ 2 , ν(dx)] heeft. Tot slot bespreken we de Esscher-transformatie die we kunnen gebruiken om de kansmaat P te transformeren naar een risiconeutrale maat Q. Veronderstel dat X = {Xt , t ≥ 0} een L´evy-proces is. De Esscher-transformatie transformeert de kansdichtheidsfunctie ft (x) van de verdeling van Xt onder P naar een nieuwe kansdichtheid ft (x, θ) met parameter θ. Definitie 17. Zij fnt (x) de kansdichtheidsfunctie, danowordt de Esscher-transformatie R∞ voor een re¨ele θ ∈ θ ∈ R| −∞ exp(θy)ft (y)dy < ∞ gedefinieerd door de nieuwe dichtheidsfunctie exp(θx)ft (x) . exp(θy)ft (y)dy −∞ ft (x, θ) = R ∞ Wanneer we de Esscher-transformatie gebruiken om te veranderen van maat, kunnen we het Radon-Nikod´ ym afgeleideproces dat correspondeert met deze verandering van maat als volgt noteren (cf. Cont and Tankov (2004)) Zt = dQ|Ft eθXt = = exp(θXt + γ(θ)t), dP|Ft E[eθXt ] (3.13) waarbij γ(θ) = − ln (E[exp(θX1 )]) de cumulantgenererende functie van X1 is (op het teken na). 3.4.4.2.1 Stochastische differentiaalvergelijking voor de temperatuur We zullen nu, zoals in het geval van de Brownse beweging, prijsformules proberen opstellen wanneer de temperatuur voorgesteld wordt via het volgende model dT (t) = dS(t) + κ(t)(T (t) − S(t))dt + σ(t)dL(t), (3.14) waarbij L een L´evy-proces is. Wanneer we X(t) = T (t) − S(t) noteren, kunnen we zoals in Sectie 2.3.2 voorgaande stochastische differentiaalvergelijking als volgt herschrijven dX(t) = κ(t)X(t)dt + σ(t)dL(t) Rt Rt ⇔d e− 0 κ(z)dz X(t) = e− 0 κ(z)dz σ(t)dL(t). 88 Hoofdstuk 3. Temperatuurderivaten prijzen Stellen we f (t, x) = e− Rt 0 κ(z)dz x, dan kunnen we de Itˆo-formule voor semimartingalen12 (zie Cont and Tankov (2004)) gebruiken om een expliciete oplossing voor de stochastische differentiaalvergelijking te bekomen. Volgens de Itˆo-formule geldt er dan f (t, X(t)) − f (0, X(0)) Z t Z t Z ∂f ∂f 1 t ∂ 2f = (u, X(u)) du + (u, X(u−)) dX(u) + (u, X(u−)) d[X, X]cu 2 ∂u ∂x 2 ∂x 0 0 0 X ∂f (u, X(u−)) . + f (u, X(u)) − f (u, X(u−)) − ∆X(u) ∂x 0≤u≤t ∆X(u) 6= 0 (3.15) Aangezien Ru Ru ∂f ∂f (u, X(u)) = −κ(u)e− 0 κ(z)dz X(u), (u, X(u−)) = e− 0 κ(z)dz en ∂u ∂x ∂ 2f (u, X(u−)) = 0 geldt, kunnen we (3.15) als volgt invullen ∂x2 e− Rt 0 κ(z)dz Z X(t) − X(0) t − Ru t Z Ru κ(u)e X(u)du + e− 0 κ(z)dz dX(u) 0 0 i X h Ru R Ru − 0 κ(z)dz − 0u κ(z)dz e + X(u) − e X(u−) − ∆X(u)e− 0 κ(z)dz . =− 0 κ(z)dz 0≤u≤t ∆X(u) 6= 0 De som in het rechterlid is gelijk aan nul, omdat ∆X(u) = X(u) − X(u−). Via (3.14) en X(t) = T (t) − S(t) bekomen we uiteindelijk e− Rt 0 κ(z)dz Z =− (T (t) − S(t)) − (T (0) − S(0)) t − κ(u)e Ru 0 κ(z)dz Z t (T (u) − S(u))du + 0 e− Ru 0 κ(z)dz [κ(u)(T (u) − S(u))du + σ(u)dL(u)]. 0 Een expliciete oplossing voor de stochastische differentiaalvergelijking (3.14) wordt dus gegeven door Rt T (t) = S(t) + e 0 κ(z)dz (T (0) − S(0)) + e Rt 0 κ(z)dz Z t σ(u)e− Ru 0 κ(z)dz dL(u). (3.16) 0 3.4.4.2.2 CAT- en PR-indices: futures en opties Z τ 2 Om de prijs van een CAT-future FCAT (t, τ1 , τ2 ) = EQ T (τ )dτ Ft te kunnen bepaτ1 len, moeten we een uitdrukking voor de cumulatieve temperatuur over het tijdsinterval 12 Als X een L´evy-proces is, dan is f (t, X(t)) geen L´evy-proces meer, het kan echter wel nog uitgedrukt worden met behulp van stochastische integralen en dus is het nog steeds een semimartingaal. Een semimartingaal is, volgens Karatzas and Kardaras (2007), een proces dat kan worden ontbonden in een term met eindige variatie en een term die een lokale martingaal is. 3.4. Dagelijkse simulatie 89 [τ1 , τ2 ] kennen. Deze uitdrukking leiden we af in Stelling 12. Stelling 12. De cumulatieve temperatuur over het tijdsinterval [τ1 , τ2 ] wordt gegeven door Z τ2 Z T (t)dt = τ1 τ2 Z 0 τ1 Rt τ1 τ1 Bewijs. Uit (3.16) volgt dat Z Z τ2 Z τ2 S(t)dt + T (t)dt = τ1 τ2 e 0 κ(z)dz (T (0) − S(0))dt S(t)dt + τ1 τ1 Z τ2 Z τ2 Z τ1 Z τ2 Rt Rt κ(z)dz u σ(u)e u κ(z)dz dtdL(u). σ(u)e dtdL(u) + + τ2 e Rt 0 κ(z)dz L(u) Z τ2 Z (T (0) − S(0))dt + σ(u)e τ1 τ1 t Rt u κ(z)dz dL(u)dt. 0 We moeten enkel de laatste term van het rechterlid nog wat herschrijven Z τ1 Z τ2 Z τ2 Z t Rt Rt κ(z)dz σ(u)e u κ(z)dz dtdL(u) σ(u)e u dL(u)dt = 0 0 τ1 Z τ2τ1Z τ2 Rt + σ(u)e u κ(z)dz dtdL(u), τ1 L(u) dit volgt rechtstreeks uit Figuur 3.1. Figuur 3.1: Verandering van integratievolgorde. Om de futureprijs expliciet te bepalen, moeten we een risiconeutrale kansmaat Q definieren. Aangezien de temperatuur niet kan worden opgeslagen, hebben we te maken met 90 Hoofdstuk 3. Temperatuurderivaten prijzen een incomplete markt en is elke equivalente maat een risiconeutrale kansmaat. In Sectie 3.4.4.1 hebben we de stelling van Girsanov gebruikt om een equivalente kansmaat Q te vinden. De stelling van Girsanov is een speciaal geval van de Esscher-transformatie wanneer de verdeling een Brownse beweging is. In het geval van een proces met sprongen, zoals een L´evy-proces, wordt de Esscher-transformatie toegepast. We gebruiken dus eigenlijk een subklasse van kansmaten aan de hand van de Esscher-transformatie. We volgen Alexandridis and Zapranis (2013) en Benth and Benth (2005) in de manier waarop zij veranderen van maat, dit is niet helemaal de standaard manier waar we reeds over gesproken hebbben in vergelijking (3.13). Veronderstel dat θ(t) een re¨eelwaardige, meetbare en begrensde functie is. Beschouw dan het volgende stochastisch proces Z t Z t Z(t) = exp θ(s)dL(s) − ϕ(θ(s))ds , 0 0 waarbij ϕ(λ) de logaritme van de momentgenererende functie, dit is de cumulantgenererende functie, van L(t) is ϕ(λ) = ln E[exp(λL(1))]. We nemen aan dat het proces Z goed gedefinieerd is onder de natuurlijke exponenti¨ele integreerbaarheidsvoorwaarden op de L´evy-maat, die we ook aannemen. We introduceren nu de kansmaat Qθ die gedefinieerd wordt via de Esscher-transformatie Qθ (A) = E[IA Z(τmax )], A ∈ F, met I de karakteristieke functie, τmax een vaste tijdshorizon die alle tijdstippen van handel voor alle relevante futurecontracten bevat en F een σ-algebra. De kansmaat Qθ is equivalent met P en we zullen in het vervolg Q noteren. Zoals in Sectie 3.4.4.1 is θ de marktprijs voor risico. Nu kunnen we de prijs voor een CAT-future bepalen. Stelling 13. De futureprijs FCAT (t, τ1 , τ2 ) op tijdstip t ≤ τ1 < τ2 wordt gegeven door FCAT (t, τ1 , τ2 ) Z τ2 Z τ2 R Z t Z τ2 Rs s κ(z)dz = S(s)ds + e0 (T (0) − S(0))ds + σ(u)e u κ(z)dz dsdL(u) τ1 τ1 Z τ2 Z τ2 Z τ1 Z 0τ1 τ1 Rs Rs + σ(u)e u κ(z)dz dsϕ0 (θ(u))du − σ(u)e u κ(z)dz dsϕ0 (θ(u))du. t L(u) t L(u) Bewijs. We bewijzen eerst dat voor een re¨eelwaardige, meetbare en begrensde functie f geldt Z EQ t τ Z f (u)dL(u)Ft = t τ f (u)ϕ0 (θ(u))du, t < τ ≤ τmax < ∞. (3.17) 3.4. Dagelijkse simulatie 91 Wanneer we gebruik maken van de formule van Bayes, de eigenschap van de onafhankelijke toenames van een L´evy-proces en de definitie van de cumulantgenererende functie van L kunnen we het volgende neerschrijven Z τ EQ f (u)dL(u)Ft tZ τ Z(τ ) =E f (u)dL(u) Ft Z(t) t Z τ Z τ d (λf (u) + θ(u))dL(u) Ft E exp = exp − ϕ(θ(u))du dλ λ=0 Zt τ Zt τ d (λf (u) + θ(u))dL(u) E exp = exp − ϕ(θ(u))du dλ t t λ=0 Z τ Z τ d ϕ(θ(u))du = exp − exp ϕ(λf (u) + θ(u))du dλ t t λ=0 Z τ f (u)ϕ0 (θ(u))du. = t In de voorlaatste gelijkheid maken we gebruik van het volgende: als g : [s, t] → R een meetbare en begrensde functie is en de voorwaarde van integreerbaarheid van de L´evy-maat geldt, dan geldt Z t Z t ϕ(g(u))du . g(u)dL(u) = exp E exp s s Het bovenstaande wordt bewezen in Cont and Tankov (2004), we laten dit bewijs achterwege. We kunnen met behulp van Stelling 12 het volgende neerschrijven Z τ EQ T (s)dsFt Z τt Z τ R s e 0 κ(z)dz (T (0) − S(0))ds S(s)ds + = t Z t Z τ t Z τZ τ Rs Rs κ(z)dz κ(z)dz + EQ σ(u)e u dsdL(u) + σ(u)e u dsdL(u)Ft 0 Z = t τ t t Z S(s)ds + Z τ Z τ Rs t L(u) κ(z)dz L(u) (T (0) − S(0))ds + EQ Rs κ(z)dz u σ(u)e dsdL(u)Ft . e + EQ t τ 0 Z t Z 0 t τ σ(u)e u κ(z)dz dsdL(u)Ft Rs 92 Hoofdstuk 3. Temperatuurderivaten prijzen Gebruik makend van (3.17) krijgen we Z τ T (s)dsFt EQ t Z τ Z τ R Z tZ τ Rs s κ(z)dz = S(s)ds + e0 (T (0) − S(0))ds + σ(u)e u κ(z)dz dsdL(u) t t 0 t Z τZ τ Rs σ(u)e u κ(z)dz dsϕ0 (θ(u))du. + t L(u) Voor de prijs van de CAT-future kunnnen we nu schrijven Z τ 2 EQ T (s)dsFt τ1 Z τ2 Z τ 1 = EQ T (s)dsFt − EQ T (s)dsFt t t Z t Z τ2 Z τ2 R Z τ2 Rs s κ(z)dz 0 σ(u)e u κ(z)dz dsdL(u) (T (0) − S(0))ds + e S(s)ds + = t t Z τ1 0 t Z τ1 R Z τ2 Z τ2 Rs s σ(u)e u κ(z)dz dsϕ0 (θ(u))du − S(s)ds − e 0 κ(z)dz (T (0) − S(0))ds + t Z tZ − Z = L(u) τ1 σ(u)e 0 τ2 Rs u κ(z)dz Z τ1 Rs t σ(u)e t τ2 t τ1 dsdL(u) − t Z Z Z tZ κ(z)dz L(u) u κ(z)dz dsϕ0 (θ(u))du L(u) S(s)ds + e0 (T (0) − S(0))ds + τ1 τ1 Z τ1 Z Z τ2 Z τ2 Rs κ(z)dz 0 σ(u)e u dsϕ (θ(u))du − + t Rs t τ2 Rs σ(u)e 0 τ1 u κ(z)dz dsdL(u) τ1 Rs σ(u)e u κ(z)dz dsϕ0 (θ(u))du, L(u) dit is de gezochte uitdrukking. Uit de definities blijkt dat de Pacific Rim index het gemiddelde van de CAT-index over een specifieke periode is. De PR-futureprijs op tijdstip t ≤ τ1 ≤ τ2 wordt gegeven door Z τ2 1 FPR (t, τ1 , τ2 ) = EQ T (s)dsFt . τ 2 − τ 1 τ1 Hierdoor vinden we dat de volgende relatie tussen een contract geschreven op de PRindex en een contract geschreven op de CAT-index geldt FPR (t, τ1 , τ2 ) = 1 FCAT (t, τ1 , τ2 ). τ2 − τ1 Wanneer we de optieprijzen expliciet willen bepalen, moeten we de volgende uitdrukking uitwerken C(t, τ, τ1 , τ2 ) = e−r(τ −t) EQ [max{F (τ, τ1 , τ2 ) − K, 0}|Ft ]. (3.18) 3.4. Dagelijkse simulatie 93 Aangezien we gebruik gemaakt hebben van een L´evy-proces om de temperatuur te modelleren, zal deze uitdrukking nog afhangen van θ(t). In dit geval introduceert het L´evy-proces een incompleetheid van de markt, die verhindert dat we de optie kunnen hedgen. Er zijn dus vele prijzen die geparametriseerd worden door θ(t). Eens θ(t) bepaald is, kan de verwachtingswaarde in (3.18) berekend worden via numerieke integratie of Monte Carlo-simulatie. Om dit te kunnen verwezenlijken, moeten we de eigenschappen van de verdeling van de stochastische variabele T (t) kennen. We kunnen de Fourier-transformatie toepassen om de dichtheid fT (x) te bepalen Z ∞ 1 fT (x) = e−isx φT (s)ds, 2π −∞ (3.19) met φT (λ) de karakteristieke functie van de stochastische variabele T (t). Wanneer de karakteristieke functie van het L´evy-proces gekend is, kunnen we de optieprijzen bepalen. Stelling 14. De karakteristieke functie van T (t) wordt onder de risiconeutrale maat Q gegeven door s≤t φT (λ) = EQ [exp(iλT (t))|Fs ] = exp(φ(λ)), waarbij Rt Z κ(z)dz φ(λ) =iλS(t) + iλe s (T (s) − S(s)) − Z t Rt ϕ iλσ(u)e u κ(z)dz + θ(u) du. + t ϕ(θ(u))du s s Bewijs. Wanneer we (3.16) substitueren in de verwachtingswaarde bekomen we voor s≤t EQ [exp(iλT (t))|Fs ] Z t Rt Rt κ(z)dz κ(z)dz dL(u) Fs = EQ exp iλ S(t) + e s (T (s) − S(s)) + σ(u)e u s Z t Rt Rt κ(z)dz κ(z)dz u s = exp iλS(t) + iλe (T (s) − S(s)) EQ exp iλ σ(u)e dL(u) Fs . s We werken de verwachtingswaarde verder uit Z t Rt κ(z)dz EQ exp iλ σ(u)e u dL(u) Fs sZ t Rt Z(t) κ(z)dz = E exp iλ σ(u)e u dL(u) Fs Z(s) s Z t Z t Z t Rt κ(z)dz u = E exp iλ σ(u)e dL(u) + θ(u)dL(u) − ϕ(θ(u))du Fs s s s 94 Hoofdstuk 3. Temperatuurderivaten prijzen Z t Z t Z t Rt κ(z)dz = exp − ϕ(θ(u))du E exp iλ σ(u)e u dL(u) + θ(u)dL(u) Fs s Z t s Zs t Rt κ(z)dz ϕ(θ(u))du E exp iλσ(u)e u + θ(u) dL(u) Fs = exp − Zs t Zs t Rt κ(z)dz u iλσ(u)e + θ(u) dL(u) ϕ(θ(u))du E exp = exp − s s Z t Z t Rt κ(z)dz ϕ(θ(u))du exp ϕ iλσ(u)e u + θ(u) du , = exp − s s de voorlaatste gelijkheid geldt door de onafhankelijke toenames van het L´evy-proces. En dus geldt er voor φ(λ) Rt κ(z)dz Z φ(λ) =iλS(t) + iλe (T (s) − S(s)) − Z t Rt κ(z)dz u + θ(u) du, ϕ iλσ(u)e + s t ϕ(θ(u))du s s waarbij ϕ de cumulantgenererende functie van L(1) is en i2 = −1 geldt. Het komt er dus op neer dat we de kansdichtheidsfunctie van het temperatuurproces moeten kennen om de optieprijzen numeriek te berekenen. Deze dichtheidsfunctie kunnen we via de Fourier-transformatie (3.19) bekomen uit de karakteristieke functie van T (t) die beschreven wordt in Stelling 14. In het geval van L´evy-processen is de momentgenererende functie (en daardoor de cumulantgenererende functie) gekend en hierdoor is ook de karakteristieke functie φ(λ), die nodig is in Stelling 14, gekend. We kennen dan uiteindelijk ook de eigenschappen van de verdeling van het gebruikte model, we kunnen immers via (3.19) de dichtheidsfunctie bepalen. Aanvankelijk hebben we dus de cumulantgenererende functie van het gebruikte L´evyproces nodig. Er bestaan veel soorten L´evy-processen, bijvoorbeeld het Poisson-proces, het Gamma-proces, het CGMY-proces, het Meixner-proces,. . . 13 Wij zijn echter ge¨ınteresseerd in de processen die we kunnen gebruiken om de temperatuur te modelleren. In Alexandridis and Zapranis (2013) werden de volgende drie processen getest om de temperatuur te beschrijven: • het Normaal Invers Gaussiaans (NIG) proces; • het getemperd stabiel proces; • het veralgemeende hyperbolisch proces. 13 Zie Schoutens (2003) voor een uitgebreide beschrijving van deze L´evy-processen. 3.4. Dagelijkse simulatie 95 Wij beperken ons hier tot de bespreking van het veralgemeende hyperbolisch proces, aangezien Alexandridis and Zapranis (2013) en Benth and Benth (2005) argumenteren dat dit een goede verdeling is om de temperatuur mee te modelleren. Volgens Schoutens (2003) kunnen we de veralgemeende hyperbolische14 (GH) verdeling via zijn karakteristieke functie als volgt defini¨eren. Definitie 18. De veralgemeende hyperbolische verdeling GH(α,β,δ,υ) wordt gedefinieerd door zijn karakteristieke functie φGH (u; α, β, δ, υ) = υ2 K δ pα2 − (β + iu)2 υ α −β p · , α2 − (β + iu)2 2 2 Kυ δ α − β 2 2 waarbij Kυ de gewijzigde Bessel-functie is. De gewijzigde Bessel-functies van de eerste soort I±υ (z) en van de derde soort Kυ (z) zijn oplossingen van de differentiaalvergelijking z2 dw d2 w +z − z 2 + υ 2 w = 0. 2 dz dz (3.20) De functie Iυ (z) kan geschreven worden als de volgende reeks Iυ (z) = ∞ z υ X 2 k=0 2 k z 4 k!Γ(υ + k + 1) , en Kυ (z) voldoet aan Kυ (z) = π Iυ (z) − I−υ (z) · , 2 sin(υπ) waarbij het rechterlid vervangen wordt door zijn limietwaarde als υ een geheel getal of nul is. De Bessel-functie Kυ kan ook in integraalvorm geschreven worden Z 1 ∞ υ−1 1 −1 Kυ (z) = u exp − z(u + u ) du, x > 0. 2 0 2 De veralgemeende hyperbolische verdeling is oneindig deelbaar en dus kunnen we een ver(GH) algemeend hyperbolisch L´evy-proces X (GH) defini¨eren. We defini¨eren X (GH) = {Xt ,t ≥ 0} als het stationair proces dat start in nul, onafhankelijke toenames heeft en waarbij (GH) de verdeling van Xt 14 de volgende karakteristieke functie heeft h i (GH) E exp iuXt = (φGH (u; α, β, δ, υ))t . In het Engels generalized hyperbolic. 96 Hoofdstuk 3. Temperatuurderivaten prijzen De L´evy-maat ν(dx) voor het veralgemeend hyperbolisch proces is nogal ingewikkeld p Z ∞ 2 exp −|x| 2y + α exp(βx) √ √ dy + υ exp(−α|x|) als υ ≥ 0, 2 y J 2 δ 2y + N 2 δ 2y |x| π 0 υ υ ν(dx) = . p Z 2 ∞ 2y + α exp −|x| exp(βx) √ √ dy als υ < 0, 2 2 |x| δ 2y + N−υ δ 2y π 2 y J−υ 0 met Jυ en Nυ de Bessel-functies uit de volgende definitie. Definitie 19. Bessel-functies van de eerste soort J±υ (z) en van de tweede soort Nυ (z) zijn oplossingen van de volgende differentiaalvergelijking z2 d2 w dw + z 2 − υ 2 w = 0. +z 2 dz dz De functie Jυ (z) kan geschreven worden als de volgende reeks 2 k −z ∞ z υ X 4 Jυ (z) = , 2 k=0 k!Γ(υ + k + 1) en Nυ (z) voldoet aan Nυ (z) = Jυ (z) cos(υπ) − J−υ (z) , sin(υπ) waarbij het rechterlid vervangen wordt door zijn limietwaarde als υ een geheel getal of nul is. Tot slot merken we nog op dat het hyperbolisch proces een speciaal geval is van het veralgemeend hyperbolisch proces. Voor υ = 1 verkrijgen we het hyperbolisch proces (HYP) (HYP) waarbij X1 de hyperbolische verdeling HYP(α,β,δ)volgt. Deze hyperbolische verdeling heeft de volgende karakteristieke functie φHYP (u; α, β, δ) = 3.4.4.2.3 21 K δ pα2 − (β + iu)2 1 α −β p · . 2 α − (β + iu)2 K1 δ α 2 − β 2 2 2 HDD- en CDD-indices: futures en opties Zoals op het einde van voorgaande sectie voor opties werd geargumenteerd, is het ook niet mogelijk om een gesloten oplossing te vinden voor futurecontracten die geschreven werden op HDD- en CDD-indices. Via Stelling 14 is het mogelijk om de karakteristieke functie van de temperatuurvariabele te bepalen, gebruik makend van de momentgenererende functie van het L´evy-proces. Via de Fourier-transformatie kunnen we dan ook de 3.4. Dagelijkse simulatie 97 dichtheidsfunctie fT (x) van de temperatuur bepalen. Deze dichtheidsfunctie kunnen we gebruiken om de prijsformule voor een CDD-future als volgt te herschrijven Z τ 2 FCDD (t, τ1 , τ2 ) = EQ max{T (s) − c, 0}dsFt Z τ2 τ1 EQ [max{T (s) − c, 0}|Ft ] ds = τ1 Z τ2 Z ∞ (x − c)fT (x)dxds. = τ1 c Op analoge manier bekomen we dat de HDD-futureprijs gegeven wordt door Z τ2 Z c (c − x)fT (x)dxds. FHDD (t, τ1 , τ2 ) = τ1 −∞ Om de optieprijzen te bepalen verwijzen we naar uitdrukking (3.18) en de bijhorende bespreking. Dit waren de prijsformules voor futures en callopties geschreven op de verschillende indices. We hebben hier geen rekening gehouden met de mogelijkheid om weersvoorspellingen te integreren in de prijs. Dit laatste bespreken we kort in de volgende sectie. 3.4.5 Prijzen met behulp van weersvoorspellingen Tot nog toe is gebleken dat het prijzen van weerderivaten sterk afhankelijk is van de beschikbare historische gegevens. Daarnaast hangt de prijs van een weerderivaat ook af van de verwachte veranderingen in het weer. Volgens Alexandridis and Zapranis (2013) kan de prijs van een weerderivaat daardoor verbeterd worden als er weersvoorspellingen beschikbaar zijn voor de periode van het contract. Deze voorspellingen zijn ook makkelijk in de prijsformules te integreren wanneer we gebruik maken van de methode die de derivaten prijst via dagelijkse simulatie. Er bestaan reeds verschillende technieken voor het ontwikkelen van weersvoorspellingen15 . Desondanks deze brede waaier van technieken moeten we er wel rekening mee houden dat het weer slechts voor enkele dagen nauwkeurig kan voorspeld worden. Maar zelfs voorspellingen op korte termijn kunnen nuttig zijn voor de prijs van contracten die reeds zijn begonnen of die in de komende dagen ingaan. Alaton et al. (2002) suggereren om weersvoorspellingen op korte termijn te gebruiken, 15 We gaan hier in deze masterproef niet dieper op in. In Alexandridis and Zapranis (2013) worden drie vaak voorkomende methodes kort besproken, met name de numerieke methoden voor weersvoorspellingen, de ensemble voorspellingen en de probabilistische voorspellingen. 98 Hoofdstuk 3. Temperatuurderivaten prijzen ze argumenteren dat voorspellingen voor meer dan een week op voorhand immers niet erg significant zijn. Om de weersvoorspellingen in de prijs (op een tijdstip dat voldoende dicht bij de start van de contractperiode ligt) te integreren, passen Alaton et al. (2002) hun model dat de temperatuur beschrijft (zie Sectie 2.4.2) aan. Ze stellen dat deze aanpassingen op verschillende manieren kunnen worden doorgevoerd. Zo kan men bijvoorbeeld, wanneer er verwacht wordt dat de temperatuur tijdens de contractperiode hoger zal zijn dan normaal, de parameter A in model (2.5) verhogen. In Ritter et al. (2011) vertrekt men van het idee dat de filtratie die gebruikt wordt in de risiconeutrale prijsformule niet alle beschikbare informatie bevat. Het is duidelijk dat deze filtratie de historische evolutie van de temperatuur bevat, maar ze bevat geen informatie over eventuele temperatuursvoorspellingen. Deze informatie is echter wel gekend door de spelers op de markt. Hiertoe stellen Ritter et al. (2011) voor om de filtratie uit te breiden tot een filtratie die alle relevante informatie bevat en dit door er de weersvoorspellingen aan toe te voegen. Daarnaast voegen ze de waarden van de voorspelde temperatuur aan de historische gegevens toe, zodat ze een model kunnen specifi¨eren dat deze voorspellingen incorporeert. De resultaten in Ritter et al. (2011) tonen aan dat de opname van weersvoorspellingen een duidelijke invloed heeft op de prijs die door de theoretische modellen gegeven wordt. Er blijkt dat, vergeleken met een referentiemodel zonder voorspellingen, de theoretische prijzen die weersvoorspellingen opnemen veel dichter bij de werkelijke futureprijzen op de markt liggen. Men merkt wel op dat dit effect slechts significant is voor de laatste twee maanden voor de vervaldatum van het contract. Het is dus zeker nuttig om voorspellingen op korte termijn te gebruiken om de prijsformules voor weerderivaten te verbeteren. Tot slot bekijken we de integratie van weersvoorspellingen nog vanuit een praktisch oogpunt. Er worden verschillende benaderingen gevolgd afhankelijk van de duur van het contract, het aantal dagen beschikbare temperatuursvoorspellingen en het aantal dagen voor de aanvang van het contract. Alexandridis and Zapranis (2013) beschrijven drie scenario’s die zich kunnen voordoen. Voor deze situaties veronderstellen we dat we beschikken over de voorspellingen van de temperatuur voor vandaag en de komende tien dagen. We hanteren ook de volgende notaties: t staat voor vandaag, τ1 voor de dag waarop het weerderivaat ingaat en τ2 voor de vervaldatum van het weerderivaat. In het concreet voorbeeld van Alexandridis and Zapranis (2013) geldt dat het contract start op 1 juli en loopt gedurende ´e´en maand. In het eerste scenario wordt de waardering van het contract uitgevoerd lang voordat het weerderivaat ingaat. Er is dus sprake van een waardering buiten de contractperiode. 3.4. Dagelijkse simulatie 99 Hierdoor zijn er op het moment van waardering nog geen weersvoorspellingen beschikbaar voor de contractperiode. We veronderstellen dat een bedrijf begin mei tot een CAT-future toetreedt. Aangezien we veronderstelden dat de beschikbare voorspellingen slechts voor de komende tien dagen zijn, dragen zij geen extra informatie bij aan de reeds bestaande modellen. De tijdspanne tussen begin mei en 1 juli is immers groter dan tien dagen. We kunnen in dit geval dus gewoon de prijsformule gebruiken die we reeds in Stelling 3 of Stelling 13 gevonden hadden. In het tweede scenario wordt de waardering van het contract opnieuw uitgevoerd vooraleer dit ingaat. Het verschil met het eerste scenario is dat de waardering geen geruime tijd voor de start van het contract plaatsvindt. Er zijn dus wel temperatuursvoorspellingen beschikbaar gedurende de contractperiode, maar niet tot helemaal aan de vervaldag. Veronderstel dus dat een bedrijf enkele dagen voor de ingang van het contract (1 juli) een CAT-future wil kopen. In dit geval kunnen we de waardering binnen de contractperiode gebruiken. Op 30 juni kent het bedrijf een voorspelling van de dagtemperatuur en van de temperatuur voor de volgende tien dagen. We kunnen de prijsformule dus verdelen in twee delen (zie Stelling 4). Het eerste deel van de formule bestaat uit de temperatuursvoorspellingen van de volgende tien dagen en is deterministisch, het is immers alsof de temperatuur voor deze tien dagen gekend is. Het tweede stuk is het stochastische deel van de contractperiode, namelijk de periode van 10 juli tot 31 juli. We kunnen de waardering van het contract als volgt noteren Z 10 FCAT (10, 1, 31) = T (s)ds + FCAT (10, 10, 31). 1 In het laatste scenario bekijken we de prijs van een CAT-future gedurende de contractperiode. Dit is dus ook het scenario waar de temperatuursvoorspellingen beschikbaar zijn tot op de vervaldag van het weerderivaat. Neem aan dat een bedrijf op 21 juli tot een CAT-future wil toetreden. Dit bedrijf kent door onze veronderstellingen de temperatuur op 21 juli en de voorspellingen voor de volgende tien dagen, d.i. tot op de vervaldatum 31 juli. Dit is eigenlijk het eenvoudigste geval. Immers, zowel de temperatuur gedurende de eerste twintig dagen van het contract als de temperatuur voor de laatste dagen is gekend. De temperatuur voor de laatste dagen van het contract is gekend omdat we deze kunnen vervangen door de voorspellingen. We kunnen opnieuw gebruik maken van de formule in Stelling 4 en verkrijgen de volgende prijsformule Z 21 FCAT (21, 1, 31) = T (s)ds + FCAT (21, 21, 31). 1 100 Hoofdstuk 3. Temperatuurderivaten prijzen De tweede term in deze gelijkheid kunnen we vervangen met behulp van de temperatuursvoorspellingen, we krijgen dan Z 21 Z FCAT (21, 1, 31) = T (s)ds + 1 31 Z T (s)ds = 21 31 T (s)ds 1 Het bovenstaande is van toepassing op de prijsformules uit Sectie 3.4.4. Er bestaan echter nog andere methoden om temperatuurderivaten te prijzen, deze alternatieve methoden maken gebruik van zogenaamde nutsfuncties. 3.5 Alternatieve methoden De markt voor weerderivaten is een incomplete markt. Hierdoor kunnen de standaard methoden die gebruikt worden om derivaten te prijzen, zoals het Black-Scholes-Mertonmodel, niet gebruikt worden. Er zijn reeds veel alternatieve methoden ontwikkeld die in een incomplete markt wel gehanteerd kunnen worden (zie Brockett et al. (2005)), bijvoorbeeld super-replicatie, kwadratische benaderingen, kwantiel hedgen en shortfall minimalisering, de aanpak via marginaal nut en indifference pricing. Sinds de modellen voor een incomplete financi¨ele markt zowel de hedgebare als de niet-hedgebare risico’s erkennen, zijn deze beter geschikt voor de waardering van weerderivaten. Sommige onderzoekers hebben dan ook geprobeerd om het probleem van de waardering van een weerderivaat te verkennen in de incomplete markt. Davis (cf. Davis (2001)) onderzoekt de aanpak via marginaal nut, gebaseerd op de veronderstelling dat agenten in de markt voor weerderivaten niet representatief zijn, maar te maken hebben met zeer specifieke weerrisico’s. Davis modelleert geaccumuleerde HDDs en grondstofprijzen via een geometrische Brownse beweging, dit leidt tot expliciete uitdrukkingen voor swap-rentes en optiewaarden. Brockett, Wang en Yang, daarentegen, passen de indifference valuation approach aan om de waardering van weerderivaten te analyseren. En Cao en Wei stellen een equilibrium valuation kader voor weerderivaten voor. Hieronder bespreken wij deze laatste twee methoden uitgebreider. 3.5.1 Indifference valuation Wij baseren deze sectie op Brockett et al. (2009). In een gemiddelde-variantie kader passen zij de indifference pricing approach aan om weerderivaten te prijzen rekening houdend met de portefeuille-effecten. De indifference pricing approach vloeit voort uit het basisbeginsel van gelijkwaardig nut. De methode is gebaseerd op argumenten van het verwachte nut en levert de 3.5. Alternatieve methoden 101 reserveringsprijzen. Ze is opgebouwd rond de voorkeuren van de beleggers in verband met risico’s die niet kunnen worden ge¨elimineerd als gevolg van de incompleetheid van de markt. De indifference pricing approach is al toegepast voor het prijzen van traditionele financi¨ele derivaten in een incomplete markt, voor het prijzen van verzekeringsproducten in de financi¨ele markt en voor het prijzen van enkele andere derivaten met exotische onderliggende. Ter illustratie van de algemene indifference pricing approach (het tijdscontinue model), veronderstellen we een dynamische markt die bestaat uit twee activa: • Een aandeel met prijsproces S = (St )0≤t≤T , met T een vaste tijdshorizon, en een spaarrekening met een constant prijsproces gelijk aan ´e´en. • Een niet-verhandelbaar actief Yt , waarop een vordering van het Europese type geschreven is. De payoff van het Europese derivaat noteren we als g(Yt ), betaalbaar op tijdstip T . Tevens wordt er geen handel in het derivaat toegestaan na inschrijving/aankoop. De individuele voorkeur voor het nemen van risico’s wordt gemodelleerd via een nutsfunctie u. In dit model wil de belegger het verwachte nut van zijn rijkdom maximaliseren door zowel rekening te houden met als zonder rekening te houden met de Europese vordering. Het optimalisatieprobleem, zonder rekening te houden met de financi¨ele vordering, is een klassiek Merton-model van optimaal investeringsgedrag, namelijk Z T V (w) = sup E u w + vt dSt , v 0 waar w staat voor de initi¨ele rijkdom van de belegger en vt voor zijn/haar investeringsstrategie (of de samenstelling van zijn/haar beleggingsportefeuille uit het risicovolle actief en het risicovrije actief op tijdstip t). V (w) is dus het maximum haalbare verwachte nut van de uiteindelijke rijkdom wanneer een belegger start met initi¨ele rijkdom w. Wanneer we rekening houden met de mogelijkheid om δ > 0 eenheden van de financi¨ele vordering te kopen/verkopen, dan wordt het optimalisatieprobleem voor de koper (b) en de verkoper (s) respectievelijk Z T b V (w) = sup E u w + vt dSt − δπ + δg(YT ) , v 0 Z T s s V (w) = sup E u w + vt dSt + δπ − δg(YT ) , b v 0 met π b en π s respectievelijk de prijzen voor het kopen en verkopen van ´e´en eenheid van de financi¨ele vordering. 102 Hoofdstuk 3. Temperatuurderivaten prijzen De verkoper zijn indifferentieprijs voor de Europese vordering g(YT ) wordt gedefinieerd door π s , deze van de koper door π b , en zodanig dat de belegger onverschillig is tegenover de volgende twee scenario’s: • Optimalisatie van zijn/haar verwachte nut zonder gebruik te maken van het derivaat. • Optimalisatie van zijn/haar verwachte nut rekening houdend met enerzijds de aansprakelijkheid (payoff) g(YT ) op vervaldatum T en anderzijds de compensatie π s (kost π b ). Daarom moeten de indifferentieprijzen π s en π b voldoen aan Z T Z T s vt dSt = sup E u w + vt dSt + δπ − δg(YT ) , sup E u w + v v 0 0 Z T Z T b sup E u w + vt dSt = sup E u w + vt dSt − δπ + δg(YT ) . v 0 v 0 Deze algemene aanpak moet nu nog toegepast worden voor weerderivaten. Door de structuur van de weerderivaten zullen hedgers deze meestal houden voor de hele looptijd van het contract (of de hedging-periode) tot op de vervaldag.16 Daarom veronderstellen Brockett et al. (2009) een twee-data-model, voor dit onderzoek is dat immers beter geschikt dan een multi-perioden-model of een tijdscontinue setting. In dit twee-data-model bestaat de financi¨ele markt uit twee activa, een risicovol actief (zoals een aandeel, de marktportefeuille,. . . ) met een stochastisch rendement r op tijdstip 1 en een spaarrekening met een bruto risicovrij rendement rf enerzijds evenals een stochastische weerindex y waarop een weerderivaat kan worden geschreven anderzijds. De payoff van dit weerderivaat wordt voorgesteld door Ry . Brockett et al. (2009) modelleren de voorkeur van de belegger voor het nemen van risico’s door de gemiddelde-variantie nutsfunctie. Deze drukt uit dat door hogere verwachte welvaart de waarde van ondernemingen stijgt, terwijl hogere volatiliteit kosten impliceert als gevolg van de toename van de kans op financi¨ele problemen en/of de effecten op toekomstige investeringsprikkels. De nutsfunctie wordt dus als volgt gedefinieerd u(x) = E(x) − λσ 2 (x), 16 Zelfs in de veronderstelling dat het weerderivaat initieel wordt gekocht/verkocht en dat de positie in dat instrument niet verandert, is het zeker mogelijk dat de belegger de gewichten van de andere effecten in de portefeuille dynamisch kan herschikken. Brockett et al. (2009) laten de invloed van de dynamische herschikking voor wat het is, zij gebruiken een statisch model omdat het wiskundig eenvoudiger en gemakkelijker te implementeren is. 3.5. Alternatieve methoden 103 met λ > 0 de risico-aversie parameter. We merken nog op dat andere doelfuncties, zoals de exponenti¨ele nutsfunctie of de macht nutsfunctie, op een analoge manier kunnen worden gebruikt. In het model van Brockett et al. (2009) wil de belegger de gemiddeldevariantie nutsfunctie van zijn/haar uiteindelijke rijkdom op tijdstip 1 maximaliseren door zowel rekening te houden met als geen rekening te houden met het weerderivaat. De koper (of hedger) wordt in het model getypeerd door een leverancier van elektriciteit met een willekeurige vraag naar stroom q en een willekeurige eenheidsprijs voor stroom p op tijdstip 1, de verkoper (of uitgever) van weerderivaten kan bijvoorbeeld een investeringsbank, een bedrijf dat energie verhandelt of een verzekeringsmaatschappij zijn. Het model incorporeert prijs risico’s (p), weer/hoeveelheid risico’s (q) en andere risico’s (r) in de financi¨ele markt. Brockett et al. (2009) waarderen weerderivaten in een hedging context waarin hedgers weerderivaten gebruiken om hun weerrisico’s te hedgen en hun nut te maximaliseren. Als eerste analyseren we de koper zijn indifferentieprijs van weerderivaten. Zonder weerhedging is het probleem van de optimale portfoliokeuze voor de koper het volgende V1b = max {u((w − a)rf + ar + pq)} , a (3.21) waarbij (w − a)rf + ar + pq de uiteindelijke rijkdom van de koper op tijdstip 1 voorstelt. Hierin is w de initi¨ele rijkdom van de belegger, a de hoeveelheid die ge¨ınvesteerd wordt in het risicovol actief, w − a de hoeveelheid die in het risicovrij actief ge¨ınvesteerd wordt en pq de inkomsten van de koper door het leveren van stroom. Rekening houdend met het gebruik van een weerderivaat om het weerrisico te hedgen, wordt het probleem van de optimale portfoliokeuze voor de koper V2b = max{u((w − δπ b )rf + a(r − rf ) + pq + δRy )}, a (3.22) met π b de koopprijs van het weerderivaat. Gebruik makend van de gemiddelde-variantie nutsfunctie u(x) = E[x] − λb σ 2 (x) kunnen we (3.21) en (3.22) expliciet bepalen. Er geldt voor (3.21) V1b = max{u((w − a)rf + ar + pq)} a = max E[(w − a)rf + ar + pq] − λb Var((w − a)rf + ar + pq) . a Wanneer we voor de eenvoud van notatie E[x] = µx , Cov(x, y) = σx,y en Var(x) = σx2 stelP P P len en aangezien E[rf ] = rf en Var( ni=1 Xi ) = ni=1 Var(Xi ) + 2 1≤i<j≤n Cov(Xi , Xj ) geldt, kunnen we dit als volgt uitschrijven V1b = max{(w − a)rf + aµr + µpq a 2 − λb ((w − a)2 σr2f + a2 σr2 + σpq + 2[(w − a)aσrf ,r + (w − a)σrf ,pq + aσr,pq ])}. 104 Hoofdstuk 3. Temperatuurderivaten prijzen Er geldt eveneens Cov(rf , x) = 0, we komen dus tot de volgende gelijkheid 2 + 2aσr,pq ) . V1b = max (w − a)rf + aµr + µpq − λb (a2 σr2 + σpq a Om het maximalisatieprobleem op te lossen, leiden we deze uitdrukking eerst af naar a en stellen we vervolgens de uitkomst gelijk aan nul, dit levert −rf + µr − λb 2aσr2 − λb 2σr,pq = 0. Wanneer we deze gelijkheid oplossen naar a bekomen we a= −rf + µr − λb 2σr,pq . 2λb σr2 Substitutie van deze waarde in de uitdrukking voor V1b levert 2 V1b = max wrf + µpq − λb σpq + a(−rf + µr − 2λb σr,pq ) − a2 λb σr2 a −rf + µr − λb 2σr,pq · (−rf + µr − 2λb σr,pq ) 2λb σr2 2 −rf + µr − λb 2σr,pq − λb σr2 2λb σr2 (−rf + µr − λb 2σr,pq )2 (−rf + µr − λb 2σr,pq )2 2 − = wrf + µpq − λb σpq + 2λb σr2 4λb σr2 (−rf + µr − λb 2σr,pq )2 2 = wrf + µpq − λb σpq + . 4λb σr2 2 = wrf + µpq − λb σpq + We kunnen nu hetzelfde doen voor (3.22) V2b = max E[(w − δπ b − a)rf + ar + pq + δRy ] − λb Var((w − δπ b − a)rf + ar + pq + δRy ) = max{(w − δπ b − a)rf + aµr + µpq + δµRy a 2 + δ 2 σR2 y + 2(aσr,pq + aδσr,Ry + δσpq,Ry ) } −λb a2 σr2 + σpq = max{(w − δπ b − a)rf + aµr + µpq + δµRy a a 2 −λb a2 σr2 − λb σpq − λb δ 2 σR2 y − 2λb aσr,pq − 2λb aδσr,Ry − 2λb δσpq,Ry }. Door het maximalisatieprobleem op te lossen, bekomen we −rf + µr − λb 2aσr2 − 2λb σr,pq − 2λb δσr,Ry = 0 −rf + µr − 2λb σr,pq − 2λb δσr,Ry ⇔ a= . 2λb σr2 3.5. Alternatieve methoden 105 Substitutie van deze waarde voor a in de uitdrukking voor V2b levert 2 V2b = max{(w − δπ b )rf + µpq + δµRy − λb σpq − λb δ 2 σR2 y − 2λb δσpq,Ry a +a(−rf + µr − 2λb σr,pq − 2λb δσr,Ry ) − a2 λb σr2 } 2 = (w − δπ b )rf + µpq + δµRy − λb σpq − λb δ 2 σR2 y − 2λb δσpq,Ry (−rf + µr − 2λb σr,pq − 2λb δσr,Ry )2 (−rf + µr − 2λb σr,pq − 2λb δσr,Ry )2 − 2λb σr2 4λb σr2 2 − λb δ 2 σR2 y − 2λb δσpq,Ry = (w − δπ b )rf + µpq + δµRy − λb σpq + (−rf + µr − 2λb σr,pq − 2λb δσr,Ry )2 + . 4λb σr2 Om de indifferentieprijs van de koper voor het weerderivaat te bepalen, stellen we V1b gelijk aan V2b , er volgt V1b = V2b 2 ⇔ wrf + µpq − λb σpq + (−rf + µr − λb 2σr,pq )2 4λb σr2 2 − λb δ 2 σR2 y − 2λb δσpq,Ry = (w − δπ b )rf + µpq + δµRy − λb σpq (−rf + µr − 2λb σr,pq − 2λb δσr,Ry )2 4λb σr2 (−rf + µr − λb 2σr,pq )2 ⇔ 4λb σr2 + = −δπ b rf + δµRy − λb δ 2 σR2 y − 2λb δσpq,Ry + (−rf + µr − 2λb σr,pq − 2λb δσr,Ry )2 . 4λb σr2 Voor de eenvoud van notatie stellen we β = −rf + µr − 2λb σr,pq , er komt dan (β − 2λb δσr,Ry )2 β2 b b 2 2 b = −δπ r + δµ − λ δ σ − 2λ δσ + pq,R f R y y Ry 4λb σr2 4λb σr2 2 β 2 − 4βλb δσr,Ry + 4(λb )2 δ 2 σr,R β2 y b b 2 2 b ⇔ b 2 = −δπ rf + δµRy − λ δ σRy − 2λ δσpq,Ry + b 2 4λ σr 4λ σr ! b 2 −βσ + λ δσ r,Ry r,Ry ⇔ 0 = δ −π b rf + µRy − λb δσR2 y − 2λb σpq,Ry + . 2 σr Deze vergelijking lossen we nu op naar π b 1 πb = rf 1 = rf µRy − λb δσR2 y − 2λb σpq,Ry + µRy − λb δσR2 y − 2λb σpq,Ry + 2 −βσr,Ry + λb δσr,R y ! σr2 2 −(−rf + µr − 2λb σr,pq )σr,Ry + λb δσr,R y σr2 ! . 106 Hoofdstuk 3. Temperatuurderivaten prijzen Gebruik makend van de gemiddelde-variantie nutsfunctie is de indifferentieprijs van de koper voor het weerderivaat gelijk aan σr,Ry (µr − rf − 2λb σr,pq − δλb σr,Ry ) 1 b b 2 b π = µRy − δλ σRy − 2λ σpq,Ry − . rf σr2 (3.23) Wegens de voorname rol van de risico-aversie parameter λ, kunnen de indifferentieprijzen van de belegger zeer subjectief zijn. In het algemeen mag de indifferentieprijs van een deelnemer niet gebruikt worden om de algemene rentabiliteit van een markt te bepalen. Echter, de transacties van weerderivaten zijn meestal over-the-counter aangepaste deals tussen een koper en een verkoper, hoewel sommige weerderivaten voor handelsdoeleinden ook worden vermeld op een aantal beurzen zoals de CME (Chicago Mercantile Exchange). De indifferentieprijzen van beleggers leveren belangrijke maatstaven aan de marktpartijen in de aangepaste markten voor weerderivaten. Om een rendabele markt te hebben, moet de marktprijs van het weerderivaat tussen de indifferentieprijs van de koper en de verkoper liggen, dit wil zeggen dat de marktprijs niet hoger is dan de indifferentieprijs van de koper. Dus, wanneer de indifference risicopremie (het verschil tussen π b en de verdisconteerde prijs aan de risicovrije rente µR y ) negatief is rf σr,Ry (µr − rf − 2λb σr,pq − δλb σr,Ry ) 1 b 2 b −δλ σRy − 2λ σpq,Ry − < 0, rf σr2 dat wil zeggen − δλb σR2 y σr2 − 2λb σpq,Ry σr2 − σr,Ry (µr − rf − 2λb σr,pq − δλb σr,Ry ) < 0 ⇔ −δλb σR2 y σr2 − 2λb σpq,Ry σr2 − σr,Ry (−2λb σr,pq − δλb σr,Ry ) < σr,Ry (µr − rf ) 2 ⇔ λb [δ(σr,R − σr2 σR2 y ) + 2(σr,Ry σr,pq − σr2 σpq,Ry )] < σr,Ry (µr − rf ), y µR y : de marktprijs van het weerderivaat (πM ) moet lager zijn dan de rf verwachte verdisconteerde17 waarde. Hieruit volgt dat de actuari¨ele prijs geen geschikte dan geldt πM < waardering van het weerderivaat is. We passen dit nu toe voor een concreet geval. Beschouw een forward contract voor het weer met tick size e 1 en uitoefenniveau K. De payoff van de koper (of hedger) voor dit contract wordt gespecifieerd door Ry = K − y, met y de weerindex. De indifferentieprijs van een forward voor het weer is dan een analogon van (3.23), waarin µRy = E[Ry ] = 17 Verdisconteren aan de risicovrije rente en onder de fysieke maat (actuari¨ele prijs). 3.5. Alternatieve methoden 107 E[K − y] = K − µy , Cov(Ry , x) = Cov(K − y, x) = −Cov(y, x) en dus Var(Ry ) = Var(y) omdat K een constante is. De prijs is dus gelijk aan 1 σr,y (µr − rf − 2λb σr,pq + δλb σr,y ) b b 2 b πf = (K − µy ) − δλ σy + 2λ σpq,y + . rf σr2 De indifference forward price F b is de waarde van K waarvoor πfb gelijk aan nul wordt, namelijk F b = µy + δλb σy2 − 2λb σpq,y − σr,y (µr − rf − 2λb σr,pq + δλb σr,y ) . σr2 Om een rendabele markt te hebben, mag de forward marktprijs niet lager zijn dan de koper zijn indifference forward price. Daarom, als de koper zijn forward indifference risicopremie (F b − µy ) positief is δλb σy2 σr,y (µr − rf − 2λb σr,pq + δλb σr,y ) > 0, − 2λ σpq,y − σr2 b dit wil zeggen δλb σy2 σr2 − 2λb σpq,y σr2 − σr,y µr − rf − 2λb σr,pq + δλb σr,y > 0 2 ⇔ λb δσy2 σr2 − 2σpq,y σr2 + σr,y 2σr,pq − δσr,y > σr,y (µr − rf ) 2 ⇔ λb δ σy2 σr2 − σr,y + 2 σr,y σr,pq − σr2 σpq,y > σr,y (µr − rf ), dan geldt FM > µy . Vandaar dat de forward marktprijs van de weerindex, FM , hoger moet zijn dan zijn verwachte waarde onder de fysieke maat. Vervolgens analyseren we de verkoper zijn indifferentieprijs van weerderivaten. Zonder rekening te houden met een weerderivaat is het probleem van de optimale portfoliokeuze voor de verkoper het volgende V1s = max{u(wrf + a(r − rf ))}. a Rekening houdend met het gebruik van een weerderivaat wordt het probleem van de optimale portfoliokeuze voor de verkoper V2s = max{u((w + δπ s )rf + a(r − rf ) − δRy )}, a met π s de verkoopprijs van het weerderivaat. We kunnen nu ook de indifferentieprijs van de verkoper berekenen gebruik makend van de gemiddelde-variantie nutsfunctie u(x) = E(x) − λs σ 2 (x). Voor V1s krijgen we V1s = max {E[wrf + a(r − rf )] − λs Var(wrf + a(r − rf ))} a 108 Hoofdstuk 3. Temperatuurderivaten prijzen = max{(w − a)rf + aµr − λs a2 σr2 } a = max{wrf + a(−rf + µr ) − λs a2 σr2 }. a Het maximalisatieprobleem oplossen naar a levert − rf + µr − 2aλs σr2 = 0 −rf + µr . ⇔a= 2λs σr2 Wanneer we deze waarde substitueren in de uitdrukking voor V1s krijgen we (−rf + µr )2 (−rf + µr )2 − 2λs σr2 4λs σr2 2 (−rf + µr ) = wrf + . 4λs σr2 V1s = wrf + Analoog geldt voor V2s V2s = max {E[(w + δπ s − a)rf + ar − δRy ] − λs Var((w + δπ s − a)rf + ar − δRy )} a = max{(w + δπ s − a)rf + aµr − δµRy − λs (a2 σr2 + δ 2 σR2 y − 2aδσr,Ry )} a = max{(w + δπ s )rf − δµRy − λs δ 2 σR2 y + a(−rf + µr + 2λs δσr,Ry ) − λs a2 σr2 }. a Wanneer we het maximalisatieprobleem oplossen, krijgen we de volgende waarde voor a − rf + µr + 2λs δσr,Ry − 2aλs σr2 = 0 −rf + µr + 2λs δσr,Ry ⇔a= . 2λs σr2 Deze waarde invullen in de uitdrukking voor V2s geeft (−rf + µr + 2λs δσr,Ry )2 (−rf + µr + 2λs δσr,Ry )2 − 2λs σr2 4λs σr2 s 2 (−rf + µr + 2λ δσr,Ry ) . + 4λs σr2 V2s = (w + δπ s )rf − δµRy − λs δ 2 σR2 y + = (w + δπ s )rf − δµRy − λs δ 2 σR2 y Om de indifferentieprijs van de verkoper voor het weerderivaat te bepalen, stellen we V1s gelijk aan V2s , er komt V1s = V2s (−rf + µr + 2λs δσr,Ry )2 (−rf + µr )2 s s 2 2 ⇔ wrf + = (w + δπ )rf − δµRy − λ δ σRy + 4λs σr2 4λs σr2 (−rf + µr )2 ⇔ 4λs σr2 3.5. Alternatieve methoden 109 = δπ s rf − δµRy − λs δ 2 σR2 y + 2 (−rf + µr )2 + 4(λs )2 δ 2 σr,R + 4(−rf + µr )λs δσr,Ry y ⇔ 0 = δπ s rf − δµRy − λs δ 2 σR2 y + s s ⇔ 0 = π rf − µRy − λ δσR2 y + 2 λs δ 2 σr,R y 2 λs δσr,R y 4λs σr2 + (−rf + µr )δσr,Ry σr2 + (−rf + µr )σr,Ry σr2 . Hieruit kunnen we π s halen 1 πs = rf µRy + λs δσR2 y − 2 λs δσr,R + (−rf + µr )σr,Ry y ! σr2 . Krachtens de gemiddelde-variantie nutsfunctie is de indifferentieprijs van de verkoper voor het weerderivaat gelijk aan σr,Ry (µr − rf + δλs σr,Ry ) 1 s s 2 π = µRy + δλ σRy − . rf σr2 (3.24) Om een rendabele markt te hebben, mag de marktprijs van het weerderivaat niet lager zijn dan deindifferentieprijs van de verkoper. Dus, als de indifference risicopremie µR π s − y van de verkoper positief is rf δλs σR2 y σr,Ry (µr − rf + δλs σr,Ry ) > 0, − σr2 dit wil zeggen δλs σR2 y σr2 − σr,Ry (µr − rf + δλs σr,Ry ) > 0 2 ⇔ λs δ(σR2 y σr2 − σr,R ) > σr,Ry (µr − rf ) y ⇔ λs > σr,Ry (µr − rf ) , 2 ) δ(σR2 y σr2 − σr,R y µR y . Op deze manier moet de marktprijs van het weerderivaat (πM ) rf hoger zijn dan de verwachte verdisconteerde18 waarde. Deze ongelijkheid geldt wanneer σr,Ry ρr,Ry < 0. Er geldt dan immers < 0 en dit impliceert σr,Ry < 0. Hieruit volgt dat σ r σ Ry dan geldt πM > σr,Ry (µr − rf ) σr,Ry (µr − rf ) = 2 2 2 σ2 δ(σRy σr − σr,Ry ) δσR2 y σr2 (1 − σ2r,Rσy2 ) Ry = 18 r σr,Ry (µr − rf ) δσR2 y σr2 (1 − ρ2r,Ry ) Opnieuw verdisconteren aan de risicovrije rente en onder de fysieke maat. 110 Hoofdstuk 3. Temperatuurderivaten prijzen negatief is. Merk op dat µr − rf positief is omdat de verwachte return van een risicovol actief r groter moet zijn dan het risicovrij rendement rf . En omdat λs > 0 verondersteld σr,Ry (µr − rf ) werd, is de ongelijkheid λs > geldig. Door de bovenstaande equiva2 δ(σR2 y σr2 − σr,R ) y µR y . Dus, wanneer de return van het risicovol actief en de lenties geldt dan ook πM > rf payoff van het weerderivaat negatief gecorreleerd zijn, is de actuari¨ele prijs geen geschikte marktprijs voor het weerderivaat. De actuari¨ele prijs is dan immers te laag. Gelijkaardige resultaten zijn van toepassing op de verkoper (of uitgever) zijn indifference forward price. De verkoper zijn indifferentieprijs van een forward voor het weer met payoff Ry = K − y kan als volgt worden uitgedrukt (analogon van (3.24)) 1 σr,y (µr − rf − λs δσr,y ) s s 2 πf = K − µy + λ δσy + . rf σr2 De indifference forward price is F s = µy − λs δσy2 − σr,y (µr − rf − λs δσr,y ) . σr2 In een rendabele markt mag de forward marktprijs niet hoger zijn dan de verkoper zijn indifference forward price. Dus, als de verkoper zijn forward indifference risicopremie (F s − µy ) negatief is −δλs σy2 σr,y (µr − rf − λs σr,y δ) − < 0, σr2 of equivalent δλs σy2 + σr,y (µr − rf − λs σr,y δ) > 0, σr2 dit wil zeggen δλs σy2 σr2 + σr,y (µr − rf − λs σr,y δ) > 0 2 ⇔ λs (δσy2 σr2 − σr,y δ) > −σr,y (µr − rf ) −σr,y (µr − rf ) 2 δ δσy2 σr2 − σr,y σr,y (µr − rf ) ⇔ λs > , 2 − σ2σ2) δ(σr,y y r ⇔ λs > dan geldt FM < µy . De forward marktprijs van de weerindex, FM , moet lager zijn dan zijn verwachte waarde onder de fysieke maat, als de markt bestaat. Deze ongelijkheid geldt wanneer ρr,y > 0. Als ρr,y > 0, dan geldt σr,y > 0. Hieruit volgt dat σr,y (µr − rf ) σr,y (µr − rf ) = 2 − σ2σ2) δ(σr,y δσy2 σr2 (ρ2r,y − 1) y r 3.5. Alternatieve methoden 111 σr,y (µr − rf ) en dus ook de onge2 − σ2σ2) δ(σr,y y r < µy geldt. Hieruit kunnen we concluderen dat wanneer de weersomstandig- negatief is. Waaruit volgt dat de ongelijkheid λs > lijkheid FM heden een positieve impact hebben op de return van het risicovol actief in de financi¨ele markt de actuari¨ele forward prijs geen geschikte forward marktprijs is in de markt voor weerderivaten. De actuari¨ele forward prijs is dan te hoog. Het is duidelijk dat een transactie in de markt voor weerderivaten enkel mogelijk is wanneer de indifferentieprijs van de koper niet lager is dan deze van de verkoper. Hierdoor is er nog een andere nodige voorwaarde voor rentabiliteit van de markt, namelijk π b ≥ π s . Deze nodige voorwaarde kan als volgt uitgedrukt worden σr,Ry (µr − rf − 2λb σr,pq − δλb σr,Ry ) 1 b 2 b µRy − δλ σRy − 2λ σpq,Ry − rf σr2 σr,Ry (µr − rf + δλs σr,Ry ) 1 µRy + δλs σR2 y − ≥ rf σr2 2 σr,R δλs σr,Ry (−2λb σr,pq − δλb σr,Ry ) y b 2 b s 2 ⇔ −δλ σRy − 2λ σpq,Ry − ≥ δλ σRy − σr2 σr2 ! ! 2 2 δσr,R σr,R δ 2σr,Ry σr,pq y y b 2 s 2 b + λ δσRy − ≥ λ δσRy − ⇔ λ −2σpq,Ry + σr2 σr2 σr2 ! 2 σr,R δ 2σr,Ry σr,pq y b s b 2 ⇔ λ −2σpq,Ry + ≥ (λ + λ ) δσRy − σr2 σr2 2 ⇔ λb −2σpq,Ry σr2 + 2σr,Ry σr,pq ≥ (λs + λb ) δσR2 y σr2 − σr,R δ y ⇔ 3.5.2 2σr,Ry σr,pq − 2σpq,Ry σr2 λb + λs . ≤ 2 λb δσR2 y σr2 − δσr,R y Equilibrium valuation Cao and Wei (2000b) veralgemenen het Lucas-model (cf. Lucas (1978)) om ook het weer te kunnen opnemen als een fundamentele variabele in de economie. Ze stellen een equilibrium valuation kader voor weerderivaten voor. Beschouw, om te beginnen, in een discrete setting een uitbreiding van de pure ruileconomie (voorgesteld door Lucas) waarin de fundamentele onzekerheden in de economie worden gedreven door twee toevalsvariabelen: • het totale dividend δ; • de weersomstandigheden Y . 112 Hoofdstuk 3. Temperatuurderivaten prijzen Het totale dividend kan gezien worden als geaggregeerde output of als dividenden op de marktportefeuille. De weersomstandigheden kunnen bepaald worden door de temperatuur, regenval, sneeuwval, vochtigheid,. . . Wij bestuderen temperatuurderivaten in het bijzonder, wat impliceert dat Y staat voor de temperatuur. De dynamiek die het totale dividend en de weervariabele bepaalt, is een exogeen proces op een gegeven kansruimte (Ω, F, P). Er is een representatieve belegger en zijn kennis zit vervat in de filtratie F = (Ft )t≥0 met Ft ≡ σ(δτ , Yτ ; τ ∈ (0, 1, 2, . . . , t)). De agent heeft een oneindige levensduur-horizon. Op de financi¨ele markt kan de representatieve agent ´e´en risicovol aandeel, discount obligaties en een eindig aantal andere voorwaardelijke vorderingen ruilen op ieder tijdstip. Het risicovol aandeel kan beschouwd worden als de marktportefeuille. Daardoor wordt de stroom van dividenden van het aandeel {δt } opgevat als het totale dividend in de economie. Het totale aanbod is genormaliseerd op ´e´en aandeel en de voorwaardelijke vorderingen zijn geschreven op het risicovol aandeel, de pure discount obligatie of de weervariabele. Het netto aanbod van alle voorwaardelijke vorderingen en de risicoloze obligatie is nul. De voorkeur van de agent wordt beschreven door een gladde tijdsadditieve verwachte nutsfunctie " V (c) = E0 ∞ X # U (ct , t) , (3.25) t=0 met U : R+ × (0, ∞) → R glad op (0, ∞) × (0, ∞) en voor elke t ∈ (0, 1, 2, . . . , ∞) is U (·, t) : R+ → R stijgend en strikt concaaf, met een continue afgeleide U 0 (·, t) op (0, ∞). Aanvankelijk wordt de agent voorzien van ´e´en risicovol aandeel. Noteer zijn posities in 0 0 de portefeuille op tijdstip t als θt = (θts , θtB , θtx ), waarbij θts , θtB en θtx het aantal aandelen is dat respectievelijk ge¨ınvesteerd is in het risicovol aandeel, de discount obligatie en andere voorwaardelijke vorderingen. Noteer de prijzen van effecten op tijdstip t door een vector Xt en de corresponderende vector van dividenden door Dt . De consumptie van de agent in de tijd wordt gefinancierd door een handelsstrategie {θt , t ≥ 0}. Zijn beslissingsprobleem bestaat er in een optimale handelsstrategie te kiezen om zo zijn verwachte nut te maximaliseren. De eerste orde voorwaarden van dit maximalisatieprobleem leveren de standaard Eulervergelijking voor consumptie op. We zullen eerst deze Eulervergelijking opstellen alvorens verder te gaan met de prijsstelling voor temperatuurderivaten die in Cao and Wei (2000b) besproken wordt. Om de Eulervergelijking voor consumptie op te stellen, volgen we de aanpak in Obstfeld 3.5. Alternatieve methoden 113 (2009) en Swisher (2012), omdat hier ook rekening wordt gehouden met de uitkering van dividenden. We gaan van start met een model dat geen rekening houdt met investeringen en/of productie. We veronderstellen eerst dat er twee periodes zijn. Op tijdstip 1 is het budget van individu i gelijk aan yi . Vanuit het standpunt op tijdstip 1 is het budget op tijdstip 2 echter een stochastische variabele. We nemen aan dat er op tijdstip 2 slechts twee mogelijke toestanden bestaan, in toestand 1 is het budget gelijk aan yi (1), in toestand 2 is het yi (2). We noteren met ci de consumptie van het individu op tijdstip 1, ci (1) en ci (2) stellen het contingentieplan voor consumptie van het individu op tijdstip 2 voor. De plannen zijn afhankelijk van de toestand die zich daadwerkelijk voordoet op tijdstip 2. De kans P dat toestand s zich voordoet is π(s), waarvoor geldt s π(s) = 1. Een belangrijke hypothese is dat het individu kiest voor het consumptieplan dat het verwachte nut U i maximaliseert. Voor U i geldt er U i = u(ci ) + β[π(1)u(ci (1)) + π(2)u(ci (2))] = u(ci ) + βE[u(ci (s))], waarbij ci (s) staat voor de consumptie in toestand s en u een nutsfunctie is. Dit is het von Neumann-Morgenstern-criterium voor het verwachte nut. Een Arrow-Debreu-effect19 voor toestand s betaalt de eigenaar ´e´en eenheid output op tijdstip 2 als toestand s optreedt en niets in alle andere gevallen20 . Zij r de intrestvoet op een obligatie, hieruit volgt dat de prijs (alle prijzen zijn in termen van consumptie op tijdstip 1) van een obligatie, die de eigenaar ´e´en eenheid output betaalt op tijdstip 2 ongeacht de toestand waarin men zich dan bevindt, gelijk is aan 1 . Vervolgens defini¨eren we de prijs op tijdstip 1 van het Arrow-Debreu-effect voor 1+r P p(s) toestand s, dit is namelijk gelijk aan . We weten eveneens dat s p(s) = 1 geldt. 1+r Stel immers dat we een Arrow-Debreu-effect voor elke toestand s zouden kopen, dan krijgen we op tijdstip 2 (ongeacht de toestand) exact ´e´en eenheid output uitbetaalt. Dit levert ons op tijdstip 2 dus hetzelfde resultaat als wanneer we een obligatie zouden bezitten, en hierdoor geldt de bovenstaande arbitrage-relatie. We bekijken het model alsof er drie goederen in bestaan, namelijk, consumptie op tijdstip 1 en consumptie op tijdstip 2 afhankelijk van de bereikte toestand. De Arrow-Debreuactivaprijzen defini¨eren de prijzen van toekomstige voorwaardelijke consumpties. Indi19 20 Ook een toestand-prijs-effect genoemd. Dit is verschillend van een risicoloze obligatie, die betaalt de eigenaar dezelfde hoeveelheid output in elke toestand. 114 Hoofdstuk 3. Temperatuurderivaten prijzen vidu i maximaliseert dus U i onderworpen aan de budgetbeperking ci + p(2) i p(1) i p(2) i p(1) i c (1) + c (2) = y i + y (1) + y (2). 1+r 1+r 1+r 1+r (3.26) Zoals in ons voorafgaand, deterministisch model spreiden mensen hun consumptie over verschillende tijdstippen (onder invloed van prijsprikkels op de verschillende tijdstippen), maar ze zijn ook van plan om hun consumptie te spreiden over verschillende toestanden afhankelijk van de prijsprikkels in de verschillende toestanden. Hoe dit in zijn werk gaat, kunnen we zien door de Lagrangiaan voor de maximalisatie van U i onderworpen aan de budgetbeperking (3.26) neer te schrijven p(1) i p(2) i i i i i i i L =U − λ c − y + c (1) − y (1) + c (2) − y (2) 1+r 1+r =u(ci ) + β(π(1)u(ci (1)) + π(2)u(ci (2))) p(1) i p(2) i i i i i i −λ c −y + c (1) − y (1) + c (2) − y (2) . 1+r 1+r De eerste orde voorwaarden worden dan verkregen door L af te leiden naar ci , ci (1), ci (2) en λi . De eerste orde voorwaarden zijn dus u0 (ci ) − λi = 0, p(1) = 0, 1+r p(2) βπ(2)u0 (ci (2)) − λi = 0, 1+r p(2) i p(1) i (c (1) − y i (1)) − (c (2) − y i (2)) = 0, − ci + y i − 1+r 1+r βπ(1)u0 (ci (1)) − λi of equivalent u0 (ci ) = λi , βπ(s)u0 (ci (s)) = λi − ci + y i = p(s) 1+r p(1) i p(2) i (c (1) − y i (1)) + (c (2) − y i (2)). 1+r 1+r Wanneer we de eerste twee voorwaarden combineren, krijgen we βπ(s)u0 (ci (s)) = u0 (ci ) p(s) , 1+r (3.27) dit is de Eulervergelijking voor een Arrow-Debreu-effect voor toestand s. Hierdoor kennen we ook de stochastische Eulervergelijking voor obligaties, we kunnen het voorgaande immers sommeren over s en bekomen dan βE[u0 (ci (s))] = u0 (ci ) 1 . 1+r 3.5. Alternatieve methoden 115 Veronderstel nu dat we over een actief beschikken dat dividenden uitbetaalt. Het actief betaalt een dividend d(s) in toestand s. Wanneer we in een model zitten dat twee perioden telt (op die manier is de waarde van het actief gelijk aan nul na de uitbetaling van het dividend), dan wordt de prijs van het actief gegeven door X p(s) d(s). 1+r s q= Gebruik makend van de Eulervergelijking voor een Arrow-Debreu-effect (3.27), kunnen we de prijs van het actief ook als volgt schrijven q= X βπ(s) s u0 (ci (s)) d(s). u0 (ci ) Wanneer we de afhankelijkheid van het individu i achterwege laten, wordt dit q=β X π(s) s u0 (c(s)) d(s), u0 (c) en dit kunnen we herschrijven tot de Eulervergelijking voor effecten qu0 (c) = β X π(s)u0 (c(s))d(s) s = βE[u0 (c(s))d(s)]. De prijs van het actief wordt dus gegeven door 0 u (c(s)) q=E β 0 d(s) . u (c) In plaats daarvan zullen we voor een actief met een lange levensduur in een economie met meer dan twee periodes hebben 0 u (ct+1 ) qt = E t β 0 (dt+1 + qt+1 ) , u (ct ) dit is een stochastische differentievergelijking in qt . Wanneer we itereren, bekomen we het volgende 0 u (ct+1 ) (dt+1 + qt+1 ) qt = Et β 0 u (ct ) 0 0 u (ct+1 ) u (ct+2 ) = Et β 0 dt+1 + Et+1 β 0 (dt+2 + qt+2 ) u (ct ) u (ct+1 ) 0 0 0 u (ct+1 ) u (ct+2 ) u (ct+2 ) = Et β 0 dt+1 + Et+1 β 0 dt+2 + Et+1 β 0 qt+2 u (ct ) u (ct+1 ) u (ct+1 ) 116 Hoofdstuk 3. Temperatuurderivaten prijzen 0 0 u0 (ct+1 ) u0 (ct+1 ) u (ct+2 ) u0 (ct+1 ) u (ct+2 ) = Et β 0 dt+1 + β 0 Et+1 β 0 dt+2 + β 0 Et+1 β 0 qt+2 u (ct ) u (ct ) u (ct+1 ) u (ct ) u (ct+1 ) 0 u0 (ct+1 ) u0 (ct+2 ) u0 (ct+1 ) u0 (ct+2 ) u (ct+1 ) dt+1 + β 0 β dt+2 + β 0 β qt+2 = Et β 0 u (ct ) u (ct ) u0 (ct+1 ) u (ct ) u0 (ct+1 ) 0 0 0 u (ct+1 ) 2 u (ct+2 ) 2 u (ct+2 ) = Et β 0 dt+1 + β dt+2 + β qt+2 u (ct ) u0 (ct ) u0 (ct ) " 2 # 0 X u0 (ct+s ) u (c ) t+2 = Et βs 0 dt+s + β 2 0 qt+2 u (c ) u (c ) t t s=1 .. . " = Et T X u0 (ct+T ) u0 (ct+s ) dt+s + β T 0 qt+T βs 0 u (ct ) u (ct ) s=1 # .. . " → Et # 0 u (c ) t+s βs 0 dt+s , u (ct ) s=1 ∞ X merk op dat wegens de toreneigenschap Et [Et+1 [Y ]] = Et [Y ] geldt en dat we gebruik gemaakt hebben van de transversale voorwaarde om de eindvoorwaarde te elimineren. Volgens Becker (2008) is de transversale voorwaarde voor een oneindig optimalisatieprobleem de randvoorwaarde die een oplossing bepaalt voor de eerste orde voorwaarden van het probleem samen met de beginvoorwaarde. De transversale voorwaarde vereist dat de actuele waarde van de toestandsvariabele naar nul convergeert als de planningshorizon naar oneindig gaat. Concreet betekent dit dus lim qt+T = 0. De eerste orde T →∞ voorwaarden samen met de transversale voorwaarden zijn voldoende om een optimum in een concaaf optimalisatieprobleem te identificeren. Volledig analoog aan het bovenstaande (met β = 1) vinden Cao and Wei (2000b) dat de prijzen van effecten op tijdstip t gegeven worden door de Eulervergelijking " ∞ # X U 0 (cτ , τ ) Xt = Et Dτ , 0 (c , t) U t τ =t+1 (3.28) met U 0 (cτ , τ ) de eerste afgeleide van de nutsfunctie naar consumptie c. Zo is de prijs van een willekeurig effect gelijk aan de som van de verwachte dividenden (Dt ), verdisconteerd aan het stochastisch marginale tarief van substitutie. In evenwicht compenseren de financi¨ele en de goederenmarkt elkaar, zodat de totale consumptie gelijk is aan de dividenden gegenereerd uit het risicovol aandeel. Er geldt dus ct = δt , waarbij ct staat voor de consumptie en δt voor de dividenden van het aandeel. Daardoor is de prijs op tijdstip t van een voorwaardelijke vordering met payoff qT op een 3.6. Conclusie 117 toekomstig tijdstip T , voorgesteld door Ct (t, T ), gelijk aan Ct (t, T ) = 1 U 0 (δt , t) Et [U 0 (δT , T )qT ], ∀t ∈ (0, T ). (3.29) Deze gelijkheid wordt verkregen uit (3.28) door slechts twee tijdstippen te beschouwen, t en T , en DT = qT te stellen. In het bijzonder is de evenwichtsprijs op tijdstip t van een risicoloze obligatie die ´e´en eenheid van consumptiegoederen betaalt op tijdstip T (dit is qT = 1 in (3.29)) en niets op alle andere tijdstippen gelijk aan B(t, T ) = 1 U 0 (δ t , t) Et [U 0 (δT , T )], ∀t ∈ (0, T ). Voorwaardelijke vorderingen op basis van een weervariabele kunnen dus gewaardeerd worden via uitdrukking (3.29) ´e´enmaal de voorkeuren van de agent, het dividendproces en het weerproces gespecifieerd zijn. Het is duidelijk dat de aanpak via maximalisering van het verwachte nut in de literatuur vaak wordt voorgesteld. Het nadeel is echter dat nutsfuncties niet voor elke belegger hetzelfde zijn. Nutsfuncties zijn immers veel te afhankelijk van de voorkeur en zijn gevoelig aan de keuze van de risico-aversie parameter. Daarnaast is het ook onmogelijk om te bepalen welke nutsfuncties er in de praktijk gebruikt worden. 3.6 Conclusie In dit hoofdstuk hebben we besproken hoe de prijs van een temperatuurderivaat bepaald kan worden. We hebben gezien dat het Black-Scholes-Merton-model in dit geval niet bruikbaar is aangezien de markt voor de temperatuur incompleet is. Deze markt is incompleet omdat de temperatuur(index) niet kan worden opgeslagen en dus geen verhandelbaar goed is. Omdat er in de praktijk toch een markt bestaat voor temperatuurderivaten, moeten er methoden bestaan om deze producten te prijzen. De drie meest gebruikte methoden zijn historical burn analysis, indexmodellering en dagelijkse simulatie. HBA prijst de derivaten aan de hand van de gemiddelde payoff van de contracten in de afgelopen n jaar. De derivaten worden dus gewaardeerd op basis van de payoff die in het verleden zou verkregen zijn. Omdat dit een zeer simpele methode is, hebben we ook eens de prijs van een fictieve calloptie berekend. Hiertoe maakten we gebruik van de Belgische temperatuurgegevens. Een nog steeds eenvoudige, maar iets nauwkeurigere aanpak is deze via indexmodellering. 118 Hoofdstuk 3. Temperatuurderivaten prijzen Hier gaat men de gebruikte temperatuurindex modelleren via parametrische verdelingen. Om de prijs van het temperatuurderivaat te bepalen, moeten de waarden van de index willekeurig uit de verdeling getrokken worden. Via deze indexwaarden kunnen vervolgens de payoffs berekend worden. Om de prijs van het derivaat te kennen, moet het gemiddelde van de payoffs genomen worden en dit resultaat moet men vervolgens verdisconteren. De meest nauwkeurige methode om temperatuurderivaten te prijzen is deze via dagelijkse simulatie. In dit geval gaat er een model opgesteld worden dat de temperatuur beschrijft op een dagelijkse basis. Zoals we in Hoofdstuk 2 reeds gezegd hebben, zijn er twee methoden om de gemiddelde dagtemperatuur te modelleren. Dit kan via een discreet (econometrische modellen) of via een continu (stochastische differentiaalvergelijking) proces. Op basis van ´e´en van de econometrische modellen die we in Hoofdstuk 2 voor de Belgische temperatuur hebben opgesteld, hebben we aangetoond hoe een optie kan geprijsd worden via discrete dagelijkse simulatie. In het geval van continue dagelijkse simulatie veronderstellen we dat de temperatuur kan gemodelleerd worden via een mean-reverting Hull-White stochastische differentiaalvergelijking. We hebben eerst aangenomen dat deze differentiaalvergelijking gedreven wordt door een Brownse beweging en nadien hebben we verondersteld dat ze gedreven wordt door een L´evy-proces. Aan de hand van deze stochastische differentiaalvergelijkingen hebben we dan wiskundige formules om de futureprijs en de prijs van een calloptie te bepalen, opgesteld. Nadien hebben we kort besproken hoe we de prijsformules nog wat scherper kunnen stellen. Dagelijkse simulatie laat immers toe om weersvoorspellingen in de prijsformules te integreren en zo de prijs van een weerderivaat te verbeteren. Dit waren de drie meest voorkomende methoden om temperatuurderivaten te prijzen. Er zijn echter ook nog vele alternatieve methoden om derivaten in een incomplete markt te prijzen. Deze methoden kunnen aangepast worden om temperatuurderivaten te prijzen, aangezien de markt voor temperatuurderivaten een incomplete markt is. Wij hebben twee van deze methoden uitgebreider besproken, de indifference valuation approach en de equilibrium valuation method. Beiden maken gebruik van nutsfuncties en stellen het maximaliseren van het verwachte nut tot doel. Bijlage A Anderstalige samenvatting For certain organizations or individuals there is a financial risk associated with adverse or unexpected weather conditions. This means that the weather becomes an additional concern for risk managers. However, in the late nineties, the financial market has thought of a solution: weather derivatives. A derivative is a financial product that derives its value from another, underlying asset. So the payoff of weather derivatives will be linked to the weather. This thesis deals with temperature derivatives in particular, so the underlying asset of these products is the temperature, or more precisely the temperature index (HDD, CDD, CAT and PR). To price a temperature derivative we need to know the temperature dynamics, which we capture in a model. There are two approaches to model the temperature. The first one uses discrete time models and the second approach makes use of continuous time models. The discrete approach specifies an econometric model, while the continuous time models for the temperature are expressed in the form of a stochastic differential equation. In this thesis we started by giving the definition of some basic econometric models and stochastic differential equations, as to understand the temperature models that are specified in the literature. To this end we deal with the AR, MA, ARMA, OrnsteinUhlenbeck, Hull-White,. . . processes. Next we have given three examples, three models that are proposed in the literature to model the temperature. The first model, by Campbell and Diebold, is an example of a discrete, econometric model. This model specifies both conditional mean dynamics and conditional variance dynamics of daily temperature, see (2.4). The second model is the famous Alaton-model. This continuous time model uses a meanreverting Ornstein-Uhlenbeck process (2.6) to model the temperature. The Benth-model uses the same Ornstein-Uhlenbeck process for the temperature, but the seasonality and 119 120 Bijlage A. Anderstalige samenvatting variance are modelled differently. We contributed to the class of discrete time models by specifying two econometric models for the Belgian temperature, based on temperature data for the period 01/01/1981 31/12/2010. These temperature data show the presence of seasonality and a (possible) trend, which we model by two different approaches, resulting in two different models. Our first approach is to model the seasonality and trend using the difference operator ∆ which is defined as ∆Yt = Yt − Yt−1 , where Yt denotes the temperature at day t. Linear trends can be eliminated by applying ∆ once. Seasonal effects of order s can be eliminated by applying the difference operator of order s, ∆s Yt = Yt − Yt−s . This approach leads us eventually to an ARMA(1,5) model for the time series ∆∆365 Yt . The first model for the Belgian temperature is therefore given by Yt =1.8097Yt−1 − 0.8097Yt−2 + Yt−365 − 1.8097Yt−366 + 0.8097Yt−367 + ut − 0.8138ut−1 − 0.2544ut−2 + 0.0168ut−3 + 0.0184ut−4 + 0.0331ut−5 , where u denotes white noise. For daily data it is more generally accepted to model the trend and seasonality using a Fourier series. We specify a Fourier series for the Belgian temperature data, using the R2 and AIC value of the specified model. We obtain the following specification of the trend and seasonality fs = a + 8 X i=1 2πs 18πs 2πs cs,i sin i + cc,i cos i + cs,9 sin . 365 365 365 The second model for the Belgian temperature is therefore given by 8 X 2πt 2πt 18πt + cc,i cos i + cs,9 sin Yt =a + cs,i sin i 365 365 365 i=1 (A.1) + θ + α1 Yt−1 + α2 Yt−2 + ut + β1 ut−1 + β2 ut−2 , of which the values for the parameters can be found in Table 2.1. Besides temperature modelling we also discuss pricing of temperature derivatives. It is almost standard procedure to talk about the Black-Scholes-Merton model when one wants to price financial derivatives. But this famous model cannot be used to price temperature derivatives, because the weather market is incomplete. However, there exist a few well discussed approaches to price temperature derivatives. One of these techniques is historical burn analysis (HBA). HBA uses the average of the 121 payoffs of the past n years for the specified contract to specify the price of the derivative. This pricing method only uses historical data and is very simplistic. We applied this HBA technique to the Belgian temperature data to price a specific call option. Another approach to price temperature derivatives is index modelling. This method models the temperature index with parametric distributions. To price a certain derivative, one should arbitrarily draw index values from the distribution and use this value to calculate the payoffs. Finally, the average of these payoffs is taken and this value is discounted. The most accurate method is the daily simulation approach. To perform this pricing method one needs a temperature model. As already said, there exist two kinds of models, discrete and continuous time models. As an own contribution we demonstrate how derivatives can be priced using a discrete time model such as (A.1). For a continuous time model (e.g. based on a mean-reverting Hull-White stochastic differential equation), we develop the pricing formulas of various temperature derivatives. First we present the formulas when the stochastic differential equation is driven by a Brownian motion. Next we present the same formulas using a stochastic differential equation but now driven by a L´evy process. We also discuss shortly how to integrate weather forecasts in the pricing formulas obtained by the continuous model. These three approaches are commonly used, but there also exist a few alternative methods. The alternative methods are designed to price derivatives in incomplete markets. The weather market is incomplete. Hence, these methods can be adapted to price temperature derivatives. We have discussed two of these approaches, the indifference valuation approach and the equilibrium valuation method. Both methods use utility functions and their goal is to maximize the expected utility of a market participant. 122 Bijlage A. Anderstalige samenvatting Bijlage B R-code B.1 Algemene eigenschappen van de tijdreeks library(forecast) library(fma) mydata = read.table("temperatuurgegevens.txt", header=TRUE) show(mydata) summary(mydata) temp<-mydata[,7] temperatuur<-ts(temp,frequency=365, start=c(1981,1)) plot.ts(temperatuur) temperatuurzondertrendenseizoen<-diff(diff(temperatuur,lag=365)) plot.ts(temperatuurzondertrendenseizoen) Acf(temperatuurzondertrendenseizoen,plot=T, main="correlogram") pacf(temperatuurzondertrendenseizoen,plot=T, main="partieel correlogram") Acf(temperatuurzondertrendenseizoen,plot=T,lag.max=800) B.2 Differentie-operator ## SARIMA-model geeft error mymodel<-arima(temperatuur,order=c(0,1,7),seasonal=c(0,1,1)) ## probeer een ARMA-model voor temperatuurzondertrendenseizoen mymodel<-arima(temperatuurzondertrendenseizoen,order=c(0,0,7)) mymodel 123 124 Bijlage B. R-code Acf(mymodel$res,plot=T,main="correlogram residuen") Box.test(mymodel$res,lag=20,type="Ljung-Box") mymodel1<-arima(temperatuurzondertrendenseizoen,order=c(1,0,1)) mymodel1 Acf(mymodel1$res,plot=T,main="correlogram residuen") Box.test(mymodel1$res,lag=20,type="Ljung-Box") mymodel2<-arima(temperatuurzondertrendenseizoen,order=c(1,0,2)) mymodel2 Acf(mymodel2$res,plot=T,main="correlogram residuen") Box.test(mymodel2$res,lag=20,type="Ljung-Box") mymodel3<-arima(temperatuurzondertrendenseizoen,order=c(2,0,2)) mymodel3 Acf(mymodel3$res,plot=T,main="correlogram residuen") Box.test(mymodel3$res,lag=20,type="Ljung-Box") mymodel4<-arima(temperatuurzondertrendenseizoen,order=c(1,0,3)) mymodel4 Acf(mymodel4$res,plot=T,main="correlogram residuen") Box.test(mymodel4$res,lag=20,type="Ljung-Box") mymodel5<-arima(temperatuurzondertrendenseizoen,order=c(1,0,4)) mymodel5 Acf(mymodel5$res,plot=T,main="correlogram residuen") Box.test(mymodel5$res,lag=20,type="Ljung-Box") mymodel6<-arima(temperatuurzondertrendenseizoen,order=c(1,0,5)) mymodel6 Acf(mymodel6$res,plot=T,main="correlogram residuen") Box.test(mymodel6$res,lag=20,type="Ljung-Box") ## vanaf ARMA(1,6) geen significante bijdrage meer mymodel7<-arima(temperatuurzondertrendenseizoen,order=c(1,0,6)) B.3. Fourierreeks mymodel7 Acf(mymodel7$res,plot=T,main="correlogram residuen") Box.test(mymodel7$res,lag=20,type="Ljung-Box") mymodel8<-arima(temperatuurzondertrendenseizoen,order=c(1,0,7)) mymodel8 Acf(mymodel8$res,plot=T,main="correlogram residuen") Box.test(mymodel8$res,lag=20,type="Ljung-Box") mymodel9<-arima(temperatuurzondertrendenseizoen,order=c(1,0,8)) mymodel9 Acf(mymodel9$res,plot=T,main="correlogram residuen") Box.test(mymodel9$res,lag=20,type="Ljung-Box") mymodel10<-arima(temperatuurzondertrendenseizoen,order=c(1,0,9)) mymodel10 Acf(mymodel10$res,plot=T,main="correlogram residuen") B.3 Fourierreeks oneyear<-rep(1:365,30) D<-365;sp<-2*pi*oneyear/D fouriermodel1<-lm(temperatuur~oneyear) plot(1:D,fouriermodel1$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") ##schatten van signaal AIC(fouriermodel1) summary(fouriermodel1) fouriermodel2<-lm(temperatuur~oneyear+sin(sp)) plot(1:D,fouriermodel2$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel2) AIC(fouriermodel2) 125 126 Bijlage B. R-code fouriermodel3<-lm(temperatuur~oneyear+sin(sp)+cos(sp)) plot(1:D,fouriermodel3$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel3) AIC(fouriermodel3) fouriermodel4<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)) plot(1:D,fouriermodel4$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel4) AIC(fouriermodel4) fouriermodel5<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp)) plot(1:D,fouriermodel5$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel5) AIC(fouriermodel5) fouriermodel6<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)) plot(1:D,fouriermodel6$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel6) AIC(fouriermodel6) fouriermodel7<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)) plot(1:D,fouriermodel7$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel7) ## cos(3sp) toevoegen niet significant, laat weg. AIC(fouriermodel7) fouriermodel8a<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) B.3. Fourierreeks 127 +sin(3*sp)+cos(3*sp) + sin(4*sp)) plot(1:D,fouriermodel8a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel8a) AIC(fouriermodel8a) fouriermodel8aa<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp) + sin(4*sp)) plot(1:D,fouriermodel8aa$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel8aa) AIC(fouriermodel8aa) fouriermodel8ab<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp)+sin(4*sp)) plot(1:D,fouriermodel8ab$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel8ab) AIC(fouriermodel8ab) ## uitlaten van niet significante sin en cos levert niet onmiddellijk ## beter model. ## R^2 en AIC zakken licht ## als we beslissen om volgende sin en cos term mee te nemen, ## laat ook voorgaande staan. fouriermodel8b<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+sin(4*sp)) plot(1:D,fouriermodel8b$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel8b) AIC(fouriermodel8b) fouriermodel8ba<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+sin(4*sp)) plot(1:D,fouriermodel8ba$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel8ba) AIC(fouriermodel8ba) 128 Bijlage B. R-code fouriermodel9<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)) plot(1:D,fouriermodel9$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel9) AIC(fouriermodel9) fouriermodel9a<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)) plot(1:D,fouriermodel9a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel9a) AIC(fouriermodel9a) fouriermodel10<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)) plot(1:D,fouriermodel10$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel10) AIC(fouriermodel10) fouriermodel10a<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+sin(4*sp)+sin(5*sp)) plot(1:D,fouriermodel10a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel10a) AIC(fouriermodel10a) fouriermodel10b<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+sin(4*sp)+sin(5*sp)) plot(1:D,fouriermodel10b$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel10b) AIC(fouriermodel10b) fouriermodel11<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)) plot(1:D,fouriermodel11$fit[1:D],type="l",xlab="s",ylab="verwachte B.3. Fourierreeks 129 temperatuur",main="Patroon Seizoenen") summary(fouriermodel11) AIC(fouriermodel11) fouriermodel11a<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)) plot(1:D,fouriermodel11a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel11a) AIC(fouriermodel11a) fouriermodel12<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp)) plot(1:D,fouriermodel12$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel12) AIC(fouriermodel12) fouriermodel13<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)) plot(1:D,fouriermodel13$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel13) AIC(fouriermodel13) fouriermodel14<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)) plot(1:D,fouriermodel14$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel14) AIC(fouriermodel14) fouriermodel15<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) 130 Bijlage B. R-code +cos(6*sp)+sin(7*sp)+cos(7*sp)) plot(1:D,fouriermodel15$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel15) AIC(fouriermodel15) fouriermodel16<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)) plot(1:D,fouriermodel16$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel16) AIC(fouriermodel16) fouriermodel16a<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)) plot(1:D,fouriermodel16a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel16a) AIC(fouriermodel16a) fouriermodel17<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+cos(8*sp)) plot(1:D,fouriermodel17$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel17) AIC(fouriermodel17) fouriermodel17a<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+cos(8*sp)) plot(1:D,fouriermodel17a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel17a) AIC(fouriermodel17a) B.3. Fourierreeks 131 fouriermodel18a<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+cos(8*sp)+sin(9*sp)) plot(1:D,fouriermodel18a$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel18a) AIC(fouriermodel18a) fouriermodel18aa<-lm(temperatuur~sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+cos(8*sp)+sin(9*sp)) plot(1:D,fouriermodel18aa$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel18aa) AIC(fouriermodel18aa) fouriermodel18b<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+sin(9*sp)) plot(1:D,fouriermodel18b$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel18b) AIC(fouriermodel18b) fouriermodel19<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+cos(8*sp)+sin(9*sp)+cos(9*sp) +sin(10*sp)+cos(10*sp)+sin(11*sp)) plot(1:D,fouriermodel19$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel19) AIC(fouriermodel19) fouriermodel20<-lm(temperatuur~oneyear+sin(sp)+cos(sp)+sin(2*sp)+cos(2*sp) +sin(3*sp)+cos(3*sp)+sin(4*sp)+cos(4*sp)+sin(5*sp)+cos(5*sp)+sin(6*sp) +cos(6*sp)+sin(7*sp)+cos(7*sp)+sin(8*sp)+cos(8*sp)+sin(9*sp)+cos(9*sp) 132 Bijlage B. R-code +sin(10*sp)+cos(10*sp)+sin(11*sp)+cos(11*sp)+sin(12*sp)+cos(12*sp) +sin(13*sp)+cos(13*sp)++sin(14*sp)+cos(14*sp)+sin(15*sp)+cos(15*sp)) plot(1:D,fouriermodel20$fit[1:D],type="l",xlab="s",ylab="verwachte temperatuur",main="Patroon Seizoenen") summary(fouriermodel20) AIC(fouriermodel20) ## residuen= temperatuur-fs restemp<-ts(fouriermodel18aa$res,frequency=365,start=c(1981,1)) plot.ts(restemp) Acf(restemp,plot=T, main="correlogram") Acf(restemp,plot=T, main="correlogram",lag.max=800) pacf(restemp,plot=T,main="partieel correlogram") pacf(restemp,plot=T,main="partieel correlogram",lag.max=800) ## model voor residuen mymodel<-arima(restemp,order=c(3,0,0)) mymodel Acf(mymodel$res,plot=T,main="correlogram residuen") Box.test(mymodel$res,lag=20,type="Ljung-Box") mymodel1<-arima(restemp,order=c(1,0,0)) mymodel1 Acf(mymodel1$res,plot=T,main="correlogram residuen") Box.test(mymodel1$res,lag=20,type="Ljung-Box") mymodel2<-arima(restemp,order=c(2,0,0)) mymodel2 Acf(mymodel2$res,plot=T,main="correlogram residuen") Box.test(mymodel2$res,lag=20,type="Ljung-Box") mymodel3<-arima(restemp,order=c(4,0,0)) mymodel3 Acf(mymodel3$res,plot=T,main="correlogram residuen") Box.test(mymodel3$res,lag=20,type="Ljung-Box") mymodel3a<-arima(restemp,order=c(5,0,0)) B.4. Prijzen op basis van discrete dagelijkse simulatie 133 mymodel3a Acf(mymodel3a$res,plot=T,main="correlogram residuen") Box.test(mymodel3a$res,lag=20,type="Ljung-Box") mymodel4<-arima(restemp,order=c(0,0,12)) mymodel4 Acf(mymodel4$res,plot=T,main="correlogram residuen") Box.test(mymodel4$res,lag=20,type="Ljung-Box") mymodel5<-arima(restemp,order=c(1,0,1)) mymodel5 Acf(mymodel5$res,plot=T,main="correlogram residuen") Box.test(mymodel5$res,lag=20,type="Ljung-Box") mymodel6<-arima(restemp,order=c(2,0,2)) mymodel6 Acf(mymodel6$res,plot=T,main="correlogram residuen") Box.test(mymodel6$res,lag=20,type="Ljung-Box") mymodel7<-arima(restemp,order=c(2,0,3)) mymodel7 Acf(mymodel7$res,plot=T,main="correlogram residuen") Box.test(mymodel7$res,lag=20,type="Ljung-Box") mymodel8<-arima(restemp,order=c(4,0,1)) mymodel8 Acf(mymodel8$res,plot=T,main="correlogram residuen") Box.test(mymodel8$res,lag=20,type="Ljung-Box") B.4 Prijzen op basis van discrete dagelijkse simulatie fs<-function(d){ sp<-2*pi*d/365; z <- 10.597201-2.572994*sin(sp)-7.375602*cos(sp) 134 Bijlage B. R-code +0.397405*sin(2*sp)-0.136062*cos(2*sp)-0.066995*sin(3*sp) +0.044629*cos(3*sp)+0.250385*sin(4*sp)+0.012405*cos(4*sp) +0.029808*sin(5*sp)+0.254952*cos(5*sp)+0.184267*sin(6*sp) -0.077406*cos(6*sp)+0.078660*sin(7*sp)-0.185825*cos(7*sp) -0.206164*sin(8*sp)+0.006139*cos(8*sp)-0.126467*sin(9*sp); z } y<-fs(121:273);y index<-vector() for (i in 1:1000){ s<-simulate(mymodel6, nsim=365, future=TRUE, seed=i) x<-s[121:273] v<-x+y cat<-sum(v) index[i]<-cat } index for(i in 1:1000){ if(index[i]-2500<=0){ index[i]=0 }else{ index[i]=(index[i]-2500)*20 } } index mean(index) Bijlage C Figuren Figuur C.1: Correlogram van de residuen van het gefitte MA(7)-model voor de tijdreeks ∆∆365 Yt . 135 136 Bijlage C. Figuren Figuur C.2: Correlogram van de residuen van het gefitte ARMA(1,1)-model voor de tijdreeks ∆∆365 Yt . Figuur C.3: Correlogram van de residuen van het gefitte ARMA(1,2)-model voor de tijdreeks ∆∆365 Yt . 137 Figuur C.4: Correlogram van de residuen van het gefitte ARMA(2,2)-model voor de tijdreeks ∆∆365 Yt . Figuur C.5: Correlogram van de residuen van het gefitte ARMA(1,3)-model voor de tijdreeks ∆∆365 Yt . 138 Bijlage C. Figuren Figuur C.6: Correlogram van de residuen van het gefitte ARMA(1,4)-model voor de tijdreeks ∆∆365 Yt . Figuur C.7: Correlogram van de tijdreeks restemp. 139 Figuur C.8: Partieel correlogram van de tijdreeks restemp. Figuur C.9: Correlogram van de residuen van het gefitte AR(5)-model voor de tijdreeks restemp. 140 Bijlage C. Figuren Figuur C.10: Correlogram van de residuen van het gefitte ARMA(1,1)-model voor de tijdreeks restemp. Bibliografie P. Alaton, B. Djehiche, and D. Stillberger. On modelling and pricing weather derivatives. Applied Mathematical Finance, 9(1):1–20, Februari 2002. A.K. Alexandridis and A.D. Zapranis. Weather Derivatives: Modeling and Pricing Weather-related Risk. Springer-Verlag, New York, 2013. R. A. Becker. Transversality condition. In S. N. Durlauf and L. E. Blume, editors, The New Palgrave Dictionary of Economics. Palgrave Macmillan, Basingstoke, 2008. J.S. Benth and F.E. Benth. Stochastic modelling of temperature variations with a view towards weather derivatives. Applied Mathematical Finance, 12(1):53–85, 2005. J.S. Benth and F.E. Benth. The volatility of temperature and pricing of weather derivatives. Quantitative Finance, 7(5):553–561, 2007. J.S. Benth and F.E. Benth. Weather derivatives and stochastic modelling of temperature. International Journal of Stochastic Analysis, 2011 Article ID 576791:21 pages, 2011. doi: 10.1155/2011/576791. J.S. Benth and F.E. Benth. A critical view on temperature modelling for application in weather derivatives markets. Energy Economics, 34(2):592–602, 2012. B.M. Bibby and M. Sørensen. Martingale estimation functions for discretely observed diffusion processes. Bernoulli, 1(1/2):17–39, 1995. T. Bollerslev and J.M. Wooldridge. Quasi-maximum likelihood estimation and inference in dynamic models with time-varying covariances. Econometric Reviews, 11(2):143– 172, 1992. A. Brix, S. Jewson, and C. Ziehmann. Weather derivative modelling and valuation: A statistical perspective. In Robert S. Dischel, editor, Climate Risk and the Weather Market, chapter 8, pages 127–150. Risk Books, June 2002. 141 142 Bibliografie P.L. Brockett, M. Wang, and C. Yang. Weather derivatives and weather risk management. Risk Management and Insurance Review, 8(1):127–140, 2005. P.L. Brockett, L.L. Golden, M. Wen, and C.C. Yang. Pricing weather derivatives using the indifference pricing approach. North American Actuarial Journal, 13(3):303–315, 2009. R. Caballero, S. Jewson, and A. Brix. Long memory in surface air temperature: Detection, modeling, and application to weather derivative valuation. Climate Research, 21 (2):127–140, 2002. S.D. Campbell and F.X. Diebold. Weather forecasting for weather derivatives. American Statistical Association, 100(469):6–16, 2005. M. Cao and J. Wei. Pricing the weather. Risk, 13(5):67–70, 2000a. M. Cao and J. Wei. Equilibrium valuation of weather derivatives. Working paper, University of Toronto & Queen’s University, Canada, May 2000b. R. Cont and P. Tankov. Financial modelling with jump processes. Chapman & Hall/CRC, London, 2004. C. Croux. Time series analysis. Cursusnota’s statistiek voor actuari¨ele wiskunde, Vrije Universiteit Brussel, Etterbeek, 2013. M. Davis. Pricing weather derivatives by marginal value. Quantitative Finance, 1(3): 305–308, 2001. U. Einmahl. Stochastische processen. Cursusnota’s stochastische processen, Vrije Universiteit Brussel, Etterbeek, 2012. E. Goetghebeur. Notes for analysis of continous data. Cursusnota’s data-analyse, Universiteit Gent, Gent, 2011. D.N. Gujarati and D.C. Porter. Basic Econometrics. McGraw-Hill, Irwin, fifth edition, 2009. C. Harris. The valuation of weather derivatives using partial differential equations. Master’s thesis, University of Reading, Reading, September 2003. S. Jewson and M. Zervos. The Black-Scholes equation for weather derivatives, August 2003. URL http://ssrn.com/abstract=436282. Bibliografie 143 R. Kaas, M. Goovaerts, J. Dhaene, and M. Denuit. Modern Actuarial Risk Theory: Using R. SpringerVerlag, Berlin Heidelberg, second edition, 2009. I. Karatzas and C. Kardaras. The num´eraire portfolio in semimartingale financial models. Finance and Stochastics, 11(4):447–493, 2007. R.E. Lucas. Asset prices in an exchange economy. Econometrica, 46(6):1429–1445, 1978. G. Meissner and J. Burke. Can we use the Black-Scholes-Merton model to value temperature options? International Journal of Financial Markets and Derivatives, 2(4): 298–313, 2011. M. Obstfeld. Optimal consumption in a frictionless world: Complete markets. Cursusnota’s Economics 202A 7, University of California, Berkeley, 2009. M. Ritter, O. Musshoff, and M. Odening. Meteorological forecasts and the pricing of temperature futures. The Journal of Derivatives, 19(2):45–60, 2011. F. Schiller, G. Seidler, and M. Wimmer. Temperature models for pricing weather derivatives. Quantitative Finance, 12(3):489–500, 2012. W. Schoutens. L´evy-Processes in Finance: Pricing Financial Derivatives. John Wiley & Sons Ltd, West Sussex, England, 2003. S.E. Shreve. Stochastic Calculus for Finance II: Continuous-Time Models. SpringerVerlag, New York, 2004. W. Smith. On the simulation and estimation of the mean-reverting Ornstein-Uhlenbeck process, 2010. URL http://commoditymodels.files.wordpress.com/2010/02/ estimating-the-parameters-of-a-mean-reverting-ornstein-uhlenbeck-process1. pdf. S. Swisher. Discussion section 4 (asset pricing with complete and incomplete markets): Answer key. Cursusnota’s Econ 714: Macroeconomic theory II, University of Wisconsin-Madison, Madison, 2012. L. Valdivieso, W. Schoutens, and F. Tuerlinckx. Maximum likelihood estimation in processes of Ornstein-Uhlenbeck type. Statistical Inference for Stochastic Processes, 12(1):1–19, 2009.
© Copyright 2024 ExpyDoc