1
INLEIDING
1
Monte-Carlo simulatie voor financiële optieprijzen Studiepunten : 2
Volg stap voor stap de tekst en los de vragen op. Bedoeling is dat je op het einde van de rit een verzorgd verslag afgeeft dat een duidelijk geheel vormt. Het verslag bevat dus niet enkel de antwoorden op de vragen, maar vormt een mooie leesbare tekst. De bijhorende MATLAB-code wordt ook afgegeven. De code moet van voldoende commentaar voorzien worden.
1
Inleiding
Stel dat u eigenaar bent van een Belgisch bedrijf dat wereldwijd producten verkoopt. Een bedrijf in Amerika plaatst een grote order om over een jaar voor 1 miljoen dollar aan producten te kopen. Als u dit bedrag over een jaar ontvangt, wilt u het graag omwisselen naar euro’s. Het is echter niet duidelijk wat op dat moment de wisselkoers dollar-euro zal zijn. Met een financiële optie kunt u zich verzekeren tegen het risico dat over een jaar de dollar-euro koers ongunstig uitvalt. Een Europese put-optie1 geeft de houder het recht, maar niet de plicht, om over T jaar een bepaald goed voor de waarde K te kunnen verkopen. De parameters T en K worden vooraf vastgelegd. In ons voorbeeld is T = 1. Als goed kiezen we het standaard bedrag van 1 dollar (de aanpassing naar andere bedragen is duidelijk). Het is natuurlijk om de waarde K in de buurt van de huidige wisselkoers S0 = 0.68 te kiezen. Als concreet voorbeeld nemen we K = 0.75. Financiële opties zijn niet gratis (waarom niet?). Cruciale vraag is wat een eerlijke prijs voor een optie is. In dit project zullen we een algemene aanpak voor het prijzen van financiële opties bekijken, toegepast op ons dollar-euro voorbeeld. Deze aanpak heet Monte–Carlo simulatie. Het wordt naast de wiskunde en economie in vele wetenschappelijke gebieden (fysica, chemie, biologie, informatica, ...) toegepast. Een belangrijk onderdeel is hierbij de numerieke oplossing van stochastische differentiaalvergelijkingen.
2
Brownse beweging
We beginnen met het implementeren van de zgn. Brownse beweging, dat een belangrijke rol in het vervolg zal spelen. 1 De term “Europese” heeft geen geografische betekenis. Het wordt in de literatuur gebruikt voor opties die alleen op het einde van de looptijd kunnen worden uitgevoerd.
3
STOCHASTISCHE INTEGRALEN
2
Een stochastisch proces is een familie U(t) (t ≥ 0) van stochastische variabelen op een gezamenlijke kansruimte (Ω, A, P). Eén trekking levert een functie van [0, ∞[ naar R. Dit noemt men een pad of realisatie. Een (standaard) Brownse beweging in R is een stochastisch proces W (t) (t ≥ 0) zodanig dat • W (0) = 0 bijna zeker • W (t) − W (s) is normaal verdeeld met verwachting 0 en variantie t − s (voor 0 ≤ s ≤ t) • voor iedere 0 = t0 < t1 < . . . < tN zijn W (tj ) − W (tj−1 ) (1 ≤ j ≤ N) onderling onafhankelijk Voor numerieke doeleinden beschouwen we een gediscretiseerde Brownse beweging, waarbij W (t) geëvalueerd wordt in discrete tijdspunten t ∈ [0, T ]. Voor gegeven geheel getal N ≥ 1, zij ∆t = T /N en tj = j · ∆t voor j = 0, 1, 2, . . . , N. Noteer de benadering van W (tj ) als Wj . Zij Zj voor j = 1, 2, . . . , N onderling onafhankelijke, normaal verdeelde stochastische veranderlijken met verwachting 0 en variantie 1. Definieer W0 = 0 en √ Wj = Wj−1 + ∆t Zj (j = 1, 2, . . . , N). (1) Dan heet W0 , W1 , W2 , . . . , WN een gediscretiseerde Brownse beweging. Vraag Hoe correspondeert (1) met de definitie van een Brownse beweging? Maak een efficiënt MATLAB-programma dat, voor willekeurig gegeven T > 0 en N ≥ 1, een gediscretiseerde Brownse beweging genereert en in een grafiek uitzet tegen t0 , t1 , t2 , . . . , tN . Verbind hierbij opeenvolgende punten met lijnstukjes. Voor een efficiënte werking is het aan te bevelen het programma te vectoriseren en van het commando cumsum gebruik te maken (hoe?). Opmerking Gebruik steeds zo weinig mogelijk lussen in je MATLAB-code. Probeer zo veel mogelijk met vectoren te werken. Iedere nieuwe uitvoering van het programma zal een nieuw plaatje geven, aangezien het een stochastisch proces betreft. Eén realisatie van een gediscretiseerde Brownse beweging (met T = 1 en N = 400) is te zien in Figuur 1. Voor de experimenten verder in de opdracht is het nuttig om de uitkomsten vergelijkbaar te maken; hiervoor stellen we de beginstaat van de ’random number generator’ in. Stellen we bijvoorbeeld randn(’state’,100), dan zullen achtereenvolgende uitvoeringen van het programma dezelfde output produceren.
3
Stochastische integralen
Gegeven een functie ϕ, dan kan de integraal N X j=1
ϕ(tj−1 )(tj − tj−1 ),
RT 0
ϕ(t)dt benaderd worden door de Riemann som (2)
4
DE EULER-MARUYAMA METHODE
3
1.5 1
W
0.5 0 −0.5 −1
0
0.2
0.4
0.6
0.8
1
t Figuur 1: Realisatie van een gediscretiseerde Brownse beweging met T = 1, N = 400.
met de discrete punten tj zoals hiervoor. De integraal is gedefinieerd door het nemen van de limiet ∆t → 0. Analoog kunnen we een som van de vorm N X
ϕ(tj−1 )(Wj − Wj−1 )
(3)
j=1
RT beschouwen, als benadering van de stochastische integraal 0 ϕ(t)dW (t). Hier integreren we ϕ met betrekking tot de Brownse beweging. Zo’n integraal heet een Itô integraal.
4
De Euler-Maruyama methode
Een scalaire, autonome stochastische differentiaalvergelijking (SDV) kan geschreven worden als Z t Z t S(t) = S0 + f (S(s))ds + g(S(s))dW (s) (0 ≤ t ≤ T ), (4) 0
0
waarbij f en g gegeven scalaire reëel-waardige functies zijn en de beginwaarde S0 een gegeven reëel getal is. Deze vorm van notatie noemen we de integraalvorm. We zullen niet verder toelichten wat het precies betekent dat S(t) een oplossing is van (4). In plaats daarvan definiëren we onmiddellijk een methode voor de numerieke oplossing van (4). Dikwijls wordt de SDV (4) herschreven in de (symbolische) differentiaalvorm dS(t) = f (S(t))dt + g(S(t))dW (t)
(0 ≤ t ≤ T ),
S(0) = S0 .
(5)
Dit is niet meer dan een compacte schrijfwijze, die we verder zullen gebruiken. Om een numerieke methode toe te passen op (5), discretiseren we eerst het interval. Stel ∆τ = T /L voor een positief geheel getal L, en τj = j · ∆τ . De numerieke benadering van S(τj ) zullen we noteren als Sj . De Euler-Maruyama (EM) methode geeft: Sj = Sj−1 + f (Sj−1 )∆τ + g(Sj−1 )(W (τj ) − W (τj−1 )),
j = 1, 2, . . . , L.
(6)
5
MONTE-CARLO SIMULATIE
4
0.9 0.85
S 0.8 0.75 0.7 0.65 0
0.2
0.4
τ
0.6
0.8
1
Figuur 2: Een realisatie van het koersverloop horende bij r = 0.06 en σ = 0.20.
Vraag Geef aan hoe men (6) in verband met (4) kan brengen. Voor de Brownse beweging gebruiken we steeds de stapgrootte ∆t = T /N. De stapgrootte ∆τ voor de numerieke methode kiezen we dan een geheel meervoud R ≥ 1 van ∆t : ∆τ = R∆t. Dit verzekert ons dat de verzameling punten {tj }, waarop het gediscretiseerde Brownse pad gebaseerd is, de punten {τj }, waarop de EM oplossing wordt berekend, bevat.
5
Monte-Carlo simulatie
We zullen de EM methode toepassen op de lineaire SDV dS(t) = r S(t)dt + σS(t)dW (t),
S(0) = S0 ,
(7)
waarbij r en σ zekere reële constantes zijn: r is de rente en σ de volatiliteit. Dit is een bekend model voor het verloop van wisselkoersen. S(t) duidt de koers aan op tijdstip t. We zullen (7) hier gebruiken als model voor de dollar-euro koers. Maak een efficiënt programma dat de dollar-euro koers simuleert op tijdstippen 0 = τ0 < τ1 < τ2 < . . . < τL = T met T = 1 en in een grafiek (τj , Sj ) uitzet verbonden met lijnstukjes. Gebruik hierbij stapgrootte ∆t = 2−8 om de Brownse beweging te simuleren en pas EM toe met stapgrootte ∆τ = R∆t met R = 4. Neem S0 = 0.68, r = 0.06 en σ = 0.20. Eén realisatie van het koersverloop is te zien in Figuur 2. We kunnen nu starten met het prijzen van de optie. De volgende formule geeft de eerlijke prijs van een Europese put-optie: V = exp(−r T ) · E[max(K − S(T ), 0)]
(8)
waarbij E de verwachting aanduidt. We kunnen de verwachting goed benaderen door een gemiddelde waarde, waarbij een groot aantal realisaties M van het koersverloop wordt gekozen: VM
M 1 X = exp(−r T ) · max(K − sm,L , 0) M m=1
(9)
6
STERKE CONVERGENTIE VAN DE EM METHODE
5
met sm,L de uitkomst voor SL in de m-de realisatie. Dit heet Monte-Carlo simulatie. Maak een efficiënt programma dat de benadering VM van de optieprijs V berekent voor willekeurige M. Pas het vervolgens toe met M = 103 , 104 , 105 om de dollar-euro put-optie te prijzen. Let op We willen hier M verschillende paden genereren van de dollar-euro koers. Dit betekent dus ook M verschillende paden van de Brownse beweging. Let daarom goed op waar je randn(’state’,100) zet, zodat je niet M keer dezelfde Brownse beweging genereert. Wel wil je dat wanneer tweemaal de benaderende optieprijs VM wordt berekend voor dezelfde waarde van M, je dezelfde uitkomst bekomt. Denk dus goed na waar je randn(’state’,100) gaat plaatsen!
6
Sterke convergentie van de EM methode
Van SDV (7) is de exacte oplossing gekend, nl. 1 S(t) = S0 exp (r − σ 2 )t + σW (t) . 2
(10)
Benader nu met Monte-Carlo simulatie de optieprijs aan de hand van de exacte oplossing (10) en vergelijk met de reeds berekende benaderingen van de optieprijs. Geef vervolgens de exacte oplossing (10) van de dollar-euro koers op [0, T ] samen met de benaderde EM oplossing in één figuur en bekijk het verschil. Wat gebeurt er als je de tijdstap bij de benaderende oplossing steeds kleiner neemt? Het voorgaande resultaat lijkt te wijzen op convergentie. Herinner dat S(τn ) en Sn stochastische variabelen zijn. Om over convergentie te kunnen spreken, is het belangrijk hoe het verschil wordt gemeten. Een methode heeft sterke orde van convergentie gelijk aan p als er een constante C bestaat zodat E|Sj − S(τj )| ≤ C ∆τ p ,
(11)
waarbij τj = j∆τ ∈ [0, T ] met j = 1, 2, 3, . . . en ∆τ voldoende klein. Als f en g voldoen aan bepaalde voorwaarden, kan aangetoond worden dat EM sterke orde van convergentie p = 12 heeft. In onze numerieke voorbeelden zullen we ons richten op de fout in het eindpunt t = T , zodus laten we strong e∆τ := E|SL − S(T )|
de EM fout in het eindpunt aanduiden in de sterke zin. Als (11) geldt met p = in [0, T ], dan geldt dit zeker op het eindpunt. We hebben dus 1
strong e∆τ ≤ C ∆τ 2
(12) 1 2
op elk vast punt
(13)
7
DE MILSTEIN METHODE
6
voor voldoende kleine ∆τ . strong Maak een efficiënt programma dat de fout e∆τ berekent voor de dollar-eurokoers S(t). Neem voor de benadering van de verwachtingswaarde het gemiddelde over 1000 gesimuleerde waardes SL . Gebruik ∆t = 2−9 voor de gediscretiseerde Brownse beweging en pas EM toe met vijf verschillende stapgroottes: ∆τ = 2m−1 ∆t met 1 ≤ m ≤ 5. Maak een loglog-grafiek waar de strong stapgrootte op de x-as staat en de bijhorende fout e∆τ op de y-as. Als de bovengrens (13) geldt met een benaderende gelijkheid, dan geldt strong log e∆τ ≈ log C +
1 log ∆τ . 2
(14)
Dit representeert een rechte met richtingscoëfficiënt 12 . Fit nu een veelterm van graad 1 (in de strong kleinste kwadraten zin) door de 5 benaderingen van de fout e∆τ die je hebt berekend. (Hiervoor bestaat een commando in MATLAB. Welk?) Wat is de richtingscoëfficiënt die je bekomt? Komt deze overeen met wat we verwachten?
7
De Milstein methode
We zagen dat de EM methode sterke orde van convergentie 21 heeft. Het is mogelijk deze sterke orde te verhogen door een correctie aan het stochastisch increment toe te voegen, Sj =Sj−1 + f (Sj−1 )∆τ + g(Sj−1 )(W (τj ) − W (τj−1 )) 1 + g(Sj−1 )g 0 (Sj−1 )((W (τj ) − W (τj−1 ))2 − ∆τ ), 2
j = 1, 2, . . . , L.
Dit heet de Milstein methode. Onderzoek met numerieke experimenten de sterke orde van convergentie van de Milstein methode. Welke orde bekom je?
8
Veralgemening model
We zagen reeds dat (7) een bekend model is voor het verloop van wisselkoersen. Dit model kan veralgemeend worden tot dS(t) = r (t)S(t)dt + σ(t)S(t)dW (t),
S(0) = S0 ,
(15)
waarbij de rente r (t) en de volatiliteit σ(t) nu afhankelijk zijn van de tijd 0 ≤ t ≤ T . Van deze SDV (15) is de exacte oplossing ook gekend, nl. Z t Z t 1 2 (r (s) − σ (s))ds + S(t) = S0 exp σ(s)dW (s) . 2 0 0
(16)
Bereken voor S0 = 0.68, r (t) = −0.01e −2t + 0.06 en σ(t) = σ = 0.2 het koersverloop en de optieprijs en de sterke orde van convergentie voor de Euler-Maruyama en de Milstein methode.
8
VERALGEMENING MODEL
7
Geef steeds de figuren of resultaten die je bekomt in je verslag weer. Bekijk ook eens de rente r (t). Is dit een logisch verloop of vind je van niet? Opmerking Om de eerlijke prijs van de Europese put-optie te bepalen, gebruiken we nu volgende formule: Z T V = exp(− r (t)dt) · E[max(K − S(T ), 0)]. (17) 0
