Technische Universiteit Delft Faculteit Elektrotechniek, Wiskunde en Informatica Delft Institute of Applied Mathematics

Technische Universiteit Delft Faculteit Elektrotechniek, Wiskunde en Informatica Delft Institute of Applied Mathematics

Eindige Differentie Methodes voor het prijzen van Europese opties met het Heston model. (Engelse titel: Finite Difference Methods for pricing European Options using the Heston model.

Verslag ten behoeve van het Delft Institute for Applied Mathematics als onderdeel ter verkrijging

van de graad van

BACHELOR OF SCIENCE in TECHNISCHE WISKUNDE

door

Pim Stuurman Delft, Nederland Augustus 2011

c 2011 door Pim Stuurman. Alle rechten voorbehouden. Copyright

BSc verslag TECHNISCHE WISKUNDE

“Eindige Differentie Methodes voor het prijzen van Europese opties met het Heston model.” (Engelse titel: “Finite Difference methods for pricing European Options using the Heston model.”)

Pim Stuurman

Technische Universiteit Delft

Begeleider Prof.dr.ir. C.W. Oosterlee

Overige commissieleden Dr.ir. F.H. van der Meulen

Dr.ir M.B. van Gijzen

Dr. G.F. Ridderbos

Augustus, 2011

Delft

Inhoudsopgave 1 Inleiding 1.1

6

Lijst van symbolen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 Introductie over optie prijs modellen

6 9

2.1

Black-Scholes model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2

Heston model voor stochastische volatiliteit ν . . . . . . . . . . . . . . . . . . . .

9

2.3

Heston-Hull-White model voor stochastische volatiliteit ν en stochastische groeifactor r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

3 Van aandeel model naar optie vergelijking 3.1

3.2

3.3

11


11

3.1.1

Afleiding van de Black-Scholes optie vergelijking . . . . . . . . . . . . . .

12

3.1.2

Exacte oplossing van Black-Scholes vergelijking . . . . . . . . . . . . . . .

14

Heston model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3.2.1

Afleiding van Heston optie vergelijking . . . . . . . . . . . . . . . . . . . .

14

Heston-Hull-White model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.3.1 3.3.2

Waarom afleiding van Heston-Hull-White optie vergelijking m.b.v. portfolio argumenten niet werkt . . . . . . . . . . . . . . . . . . . . . . . . . .

17

Heston-Hull-White optie vergelijking . . . . . . . . . . . . . . . . . . . . .

18

4 Details van het Heston model

19

4.1

Motivatie voor Heston stochastische volatiliteit model . . . . . . . . . . . . . . .

19

4.2

Keuze van g(S, ν, t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

4.3

Co¨ ordinatentransformatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

4.3.1

Co¨ ordinatentransformatie S → ln(S), ν → ln(ν) . . . . . . . . . . . . . . .

20

4.3.2

Algemene co¨ ordinaten transformatie . . . . . . . . . . . . . . . . . . . . .

22

Feller conditie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.4

5 Numerieke implementatie

25

5.1

Eindige differenties voor Black-Scholes . . . . . . . . . . . . . . . . . . . . . . . .

25

5.2

Eindige differenties voor Heston . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

4

6 Resultaten

30

6.1


30

6.2

Heston model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

7 Conclusie

44

8 Referenties

44

A Itˆ o’s lemma

45

A.1 Itˆ o’s lemma voor 1 variabele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

A.2 Itˆ o’s lemma voor 3 variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

B Put-Call pariteit

47

C Upwind schema

48

D MATLAB Appendix

48

D.1 BlackScholesCall.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

D.2 BSExactCallPlot.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

D.3 BlackScholesExactCall.m

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

D.4 TestBlackScholes.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

D.5 HestonCall.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

D.6 HestonCallIte.m

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

D.7 HestonExactCall.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

D.8 TestHeston.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

5

1

Inleiding

In de financiële wereld wordt veel gehandeld in producten die opties genoemd worden. In dit onderzoek zullen we een aantal modellen analyseren die gebruikt worden om de waarde van deze producten te bepalen. Deze modellen gaan uit van stochastische differentiaalvergelijkingen om het gedrag van een aandeel te beschrijven. In sectie 2 worden opties ge¨ıntroduceerd en laten we drie verschillende modellen voor het gedrag van aandelen zien: het Black-Scholes, Heston en Heston-Hull-White model. Vervolgens leiden we in hoofdstuk 3 uit deze modellen partiële differentiaalvergelijkingen af die de waarde van een optie bepalen. Voor het Heston-Hull-White model blijkt de aanpak die bij de andere twee modellen succesvol was, niet direct te werken. De reden hiervoor wordt toegelicht. Er blijkt wel een andere methode te zijn om uit het Heston-Hull-White model een vergelijking af te leiden. In de rest van het verslag werken we verder met de modellen van Black-Scholes en Heston. In sectie 4 worden een aantal details van het Heston model besproken. Zo tonen we aan dat er geen directe co¨ ordinatentransformatie bestaat die de numerieke implementatie vergemakkelijkt. De modellen worden numeriek ge¨ımplementeerd in hoofdstuk 5. We maken hiervoor gebruik van eindige differenties. In sectie 6 worden de resultaten voor de modellen getest. Hierbij ontdekken we dat de randvoorwaarden die aan het begin van het verslag zijn ge¨ıntroduceerd, problemen opleveren bij numerieke implementatie. Vervolgens leiden we nieuwe randvoorwaardes af, die beter presteren. Het hoofdresultaat van het onderzoek is een verzameling randvoorwaarden die goed functioneren bij het numerieke schema. In sectie 7 worden deze en andere conclusies gepresenteerd. In de Appendix zijn overige wiskundige details en MATLAB programma’s te vinden.

1.1

Lijst van symbolen

In onderstaand overzicht worden alle symbolen uit het verslag toegelicht. In veel gevallen worden de symbolen in het verslag nog verder toegelicht. • S: waarde van een aandeel. • r: groeifactor van de waarde van een aandeel (of rente). • ν: volatiliteit van een aandeel. • t: tijdstip. • W : Wiener proces (ook bekend als Brownse beweging). • ν: theoretische verwachting voor de volatiliteit ν. • κ: ’mean reversion speed’ voor volatiliteit ν (snelheid waarmee de volatileit ν terugkeert naar zijn theoretische verwachting ν). 6

• γ: volatiliteit van de volatiliteit (ook bekend als tweede-orde volatiliteit). • ρ: correlatie-coëfficient tussen twee Wiener processen. √ • σ: wortel van de volatiliteit ( ν). • V : waarde van een optie. • C: waarde van een Call-optie. • P : waarde van een Put-optie. • K: uitvoerprijs van een optie (ook bekend als strike). • Cmarkt : prijs van een Call-optie in de markt. • CBS : waarde van een Call-optie volgens Black-Scholes model. • T : expiratietijd. • χ2 (x, y): non-centrale Chi-kwadraat verdeling met graad x en non-centraliteit parameter y. • q: maat voor de Feller-conditie. Als q > 0 is aan de Feller-conditie voldaan, als −1 ≤ q ≤ 0 is niet aan de conditie voldaan. • r: theoretische verwachting voor de groei van een aandeel. • λ: ’mean reversion speed’ voor de groei van een aandeel. (snelheid waarmee de groei r van een aandeel terugkeert naar zijn theoretische verwachting r. • η: volatiliteit van de groei van een aandeel. • PDV: Partiële DifferentiaalVergelijking. • G: stochast die voldoet aan de eis van Itô’s lemma. • Π: waarde van een portfolio. • ∆: hoeveelheid (bijvoorbeeld aandelen of opties). • Smax : maximale waarde van S (groot t.o.v. K, Smax >> K). • N (x): verdelingsfunctie van een standaardnormaal verdeelde stochast (N (0, 1)) in het punt x. • ln(x): natuurlijke logaritme in het punt x. • g: functie die gebruikt wordt bij scheiding van variabelen in het Heston model. • φ: marktprijs van volatiliteitsrisico. • θ: lineaire coëfficient tussen φ en ν. • κ0 : risico-aangepaste ’mean reversion speed’ voor volatiliteit. • ν 0 : risico-aangepaste theoretische verwachting voor volatiliteit ν.

7

ˆ getransformeerde S. • S: • νˆ: getransformeerde ν. • f (S): transformatiefunctie voor S. • u(ν): transformatiefunctie voor ν. • ai , bi : parameters van stochast Si in afleiding van Itô’s lemma voor 3 stochasten. • τ : tijd tot expiratiedatum (τ = T − t). • NS : aantal punten bij discretisatie in de aandeelwaarde S. • Nτ : aantal punten bij discretisatie in de tijd tot expiratiedatum τ . • Nν : aantal punten bij discretisatie in de volatiliteit ν. • Sj : aandeelwaarde op gridpunt j. • τi : tijd tot expiratiedatum in gridpunt i • νm : volatiliteit in gridpunt m. • k: stapgrootte in de tijd tot expiratiedatum τ . • h: stapgrootte in de aandeelwaarde S. • l: stapgrootte in de volatiliteit ν. • Vji : optiewaarde op gridpunt (Sj , τi ). i : optiewaarde op gridpunt (S , ν , τ ). • Vj,m j m i

• Vi : vector met de numerieke oplossing op tijdstip τi . • F : matrix van numerieke schema (tijdstipsonafhankelijk). • pi : vector voor numerieke schema (tijdstipsafhankelijk). • I: identiteitsmatrix. • Di : submatrix van numerieke schema. • Ti : submatrix van numerieke schema. • n: indicator bij horizontale nummering in ν-richting. • C: constante bij de stabiliteitsanalyse van Black-Scholes.

8

2

Introductie over optie prijs modellen

Opties zijn financiële producten die hun eigenaar het recht (maar niet de verplichting) geven een aandeel op een afgesproken datum tegen een afgesproken prijs te kopen/verkopen. In het geval van een kopen spreken we van een Call optie, bij verkopen noemen we dit een Put optie. Wanneer alleen op de afgesproken datum de optie uitgevoerd kan worden (en dus niet tussentijds), spreken we van Europese opties. Voordat we een model voor de prijzen van opties kunnen ontwikkelen, hebben we eerst een model voor de onderliggende aandelen nodig. We beschrijven drie verschillende modellen. Een veel gebruikt model voor aandelenontwikkeling is dat van een Stochastische DifferentiaalVergelijking. Dit is een (Partiële) Differentiaal Vergelijking, waarbij één of meerdere uitdrukkingen random processen bevat.

2.1

Black-Scholes model

Het beroemdste model voor aandeelwaardes is het zogenaamde Black-Scholes model. √ dS = rdt + νdW, S

(1)

waarbij S de aandeelprijs is, r de rate of return (of rente), ν de volatiliteit en W een Wiener proces is (ook bekend als Brownse beweging). We kunnen het model in stappen opbouwen. Allereerst beschouwen we een verandering in de aandeelwaarde dS in een tijdstap dt. Voor kleine dt, nemen we aan dat de groei van de aandeelwaarde lineair afhangt van de huidige aandeelwaarde met groeifactor r: dS = rS. dt Uit de praktijk weten we dat de grafiek van een aandeelwaarde een random gedrag vertoont (zie figuur 1). Om dit fenomeen in ons model te krijgen, voegen we een ruisterm toe. Dit doen we met behulp van een zogenaamd Wiener proces. Als we hier een mate van schokkerigheid (de volatiliteit) ν aan toe voegen, komen we uit op het Black-Scholes model voor aandelen 1. √ Merk op dat ν een keuze is die we vooral gebruiken voor het Heston model dat we later zullen afleiden. In veel literatuur wordt in het Black-Scholes model gebruik gemaakt van de √ parameterσ = ν. In onze numerieke analyse van het Black-Scholes model zullen wij dit ook doen.

2.2

Heston model voor stochastische volatiliteit ν

Het Heston model is een uitbreiding op het Black-Scholes model, waarbij we de volatiliteit ν ook als een stochastische differentiaalvergelijking beschrijven. In sectie 4.1 wordt deze keuze verder gemotiveerd.

9

Figuur 1: Boven: Waarde van Microsoft aandeel in $ van mei 2010 t/m mei 2011. We zien duidelijk de randomness van de grafiek. Onder: Hoeveelheid beschikbare Microsoft aandelen in dezelfde periode. We zullen de variantie ν uit (1) beschouwen als een mean reverting square root process” [2]. ” Dit betekent dat de volatiliteit zich stochastisch zal bewegen rond een theoretische verwachte waarde. Het model wordt dan

dS = rSdt + dνt

√

νSdW1 , √ = κ(ν − ν)dt + γ νdW2 ,

(2)

waarbij κ > 0 de snelheid is waarmee de volatiliteit ν terugkeert naar zijn theoretische verwachting ν > 0 en γ > 0 de tweede-orde volatiliteit (variantie van de volatiliteit). W is wederom een Wiener proces. De Wiener processen W1 en W2 zijn gecorreleerd: dW1 dW2 = ρdt, de overige variabelen zijn zoals hierboven beschreven.

2.3

Heston-Hull-White model voor stochastische volatiliteit ν en stochastische groeifactor r

Een verdere uitbreiding van het Heston model kan gemaakt worden door ook de groeifactor r met een stochastisch proces te beschrijven. Hiervoor wordt een zogenaamd “extended OrsteinUhlenbeck mean reverting process”gebruikt, zie bijvoorbeeld [6]. Dit proces is vergelijkbaar met het proces voor de volatiliteit, maar heeft een andere factor voor de ruisterm. Het aandeelmodel wordt nu 10

dS = rSdt + dνt

√

νSdW1 , √ = κ(ν − ν)dt + γ νdW2 ,

drt = λ(r − r)dt + ηdW3 ,

(3) (4)

waarbij λ > 0 de snelheid is waarmee de interest rate r terugkeert naar zijn theoretische verwachting r > 0 en η de algemene volatiliteit bepaalt. W3 is wederom een Wiener proces en de Wienerprocessen zijn onderling gecorreleerd met dWi dWj = ρij .

3

Van aandeel model naar optie vergelijking

In deze sectie leiden we voor alle drie modellen een Partiële DifferentiaalVergelijking (PDV) af om optiewaarden te berekenen. Voor het Black-Scholes en Heston model doen we dit met behulp van portfolio argumenten en Itø’s Lemma. Vervolgens laten we zien waarom dit voor het Heston-Hull-White model niet eenvoudig mogelijk is.

3.1

Black-Scholes model

√ In het Black-Scholes model beschouwen we een aandeel-model, waarbij r en σ = ν als bekende parameters (of functies) worden gegeven. We hebben dan vergelijking (1) voor het onderliggende √ aandeel. σ is de wortel van de variantie ( ν in het oorspronkelijke model). In het Black-Scholes model is het makkelijker om deze notatie aan te houden. Voordat we de Black-Scholes afleiding presenteren, is het belangrijk om kort op te merken welke aannames noodzakelijk zijn. De volgende lijst komt uit [9]: • Er zijn geen transactiekosten. • Het onderliggende aandeel betaalt geen dividend uit gedurende de looptijd van de optie. • Er zijn geen arbitragemogelijkheden. • Het onderliggende aandeel kan continu verhandeld worden. • Short selling is toegestaan. Enkele termen uit deze lijst vereisen een korte toelichting. Een arbitragemogelijkheid treedt op als er een investering bestaat die meer oplevert dan de rente op een spaarrekening, zonder enig risico toe te voegen. Equivalent is te zeggen dat elke risicoloze investering evenveel waard moet zijn als geld op een spaarrekening zetten. Short selling betekent dat er aandelen verkocht mogen worden, die de eigenaar nog niet in bezit heeft, maar later zou kunnen kopen.

11

3.1.1

Afleiding van de Black-Scholes optie vergelijking

Bij de Black-Scholes afleiding zullen we gebruik maken van Itô’s lemma, een bekende stelling uit de stochastische analyse. De stelling wordt toegelicht in Appendix A.1. We beschouwen een optie (put of call maakt op dit moment niet uit) met een waarde V (S, t) die alleen afhankelijk is van S en t. Uit Itˆ o’s lemma voor tijdsafhankelijke functies, zie vergelijking (67), volgt dan

dV = σS

∂V ∂V ∂V 1 ∂2V dt. dW + rS + σ2S 2 2 + ∂S ∂S 2 ∂S ∂t

(5)

Nu construeren we een portfolio Π bestaande uit een optie V en een aantal −∆ onderliggende aandelen. We zullen later zien wat dit aantal is. De waarde van dit portfolio is Π = V − ∆S De verandering in waarde van dit portfolio per tijdstap is dΠ = dV − ∆dS waarbij ∆ constant gehouden wordt gedurende de tijdstap. Nu combineren we (1) en (5), zodat dΠ = σS

∂V ∂V 1 2 2 ∂2V ∂V − ∆ dW + rS + σ S + − r∆S dt. ∂S ∂S 2 ∂S 2 ∂t

(6)

∂V We zien dat de stochastische term σS ∂V ∂S − ∆ dW wegvalt als ∆ = ∂S . Merk op dat we deze waarde alleen hebben aan het begin van elke tijdstap dt. Substitueren van ∆ = ∂V ∂S in (6) levert dΠ =

1 2 2 ∂2V ∂V σ S + 2 ∂S 2 ∂t

dt.

(7)

We kunnen (7) verder behandelen met behulp van enkele arbitrageargumenten. In een tijdstap dt, zou de return voor een bedrag Π in aandelen een groei van rΠdt maken. Als de rechterkant van (7) groter dan rΠdt zou zijn, kan een investeerder risicoloos geld verdienen door Π te lenen en te investeren in het portfolio. Als de rechterkant (7) kleiner dan rΠdt is, kan de investeerder short gaan op het portfolio en Π op een spaarrekening plaatsen. Beide gevallen gaan dus tegen de arbitrageaanname in, waardoor we kunnen concluderen dat rΠdt =

1 2 2 ∂2V ∂V σ S + 2 2 ∂S ∂t

dt.

(8)

Door substitutie van (Π = V − ∆S) en links en rechts te delen door dt, vinden we de befaamde Black-Scholes formule voor opties:

12

∂V 1 ∂2V ∂V + σ 2 S 2 2 + rS − rV = 0 ∂t 2 ∂S ∂S

(9)

Om de Black-Scholes vergelijking op te lossen zijn er randvoorwaarden nodig. De formule zoals gegeven in (9) geldt voor zowel put als call opties, maar ze hebben verschillende randvoorwaarden. We zullen eerst de call opties (genoteerd als C(S, t)) behandelen. Vanwege de eerste orde afgeleide in van C in de tijd t, hebben we 1 begin- of eindvoorwaarde nodig. We kunnen een eindvoorwaarde vinden door een arbitrageargument te gebruiken [1]. We weten de waarde van een call optie op expiratietijd T voldoet aan C(S, T ) = max(S(T ) − K, 0).

(10)

Deze conditie gebruiken we als eindvoorwaarde. Omdat we een tweede orde afgeleide in de S-ruimte hebben, zoeken we twee randvoorwaarden in S. We weten dat als een aandeel S(t) op een tijdstip t gelijk is aan 0, de waarde van het aandeel 0 zal blijven. Hieruit volgt dat de call ook waarde 0 krijgt. Zo vinden we de eerste randvoorwaarde C(0, t) = 0

(11)

voor alle t. Op vergelijkbare wijze zien we dat als een aandeelprijs zeer groot wordt (S >> K), het aandeel een zeer grote waarde zal behouden, waardoor de strike K insignificant wordt. De waarde van een Call-optie wordt dan simpelweg de waarde van het aandeel. We zien dus, voor Smax >> K, C(Smax , t) ≈ Smax .

(12)

Analoog aan deze redenering vinden we voor een put optie P (S, t) de eindvoorwaarde P (S, T ) = max(K − S(T ), 0).

(13)

Voor de randvoorwaarden zien we dat voor zeer grote aandeelwaardes S >> K de put-optie nooit in the money zal komen, waardoor de optie 0 waard is. We zien dus, voor Smax >> K, P (Smax , t) ≈ 0.

(14)

Tot slot zien we dat als een aandeel op een bepaald tijdstip 0 waard is, het aandeel 0 waard blijft, waaruit volgt dat P (0, T ) = K op expiratietijdstip. We kunnen dit met de rente r verdisconteren, zodat we vinden P (0, t) = Ke−r(T −t) .

(15)

De waarde van een put-optie kan ook gerelateerd worden aan de waarde van een call-optie via de put-call pariteit. Voor meer informatie hierover, zie appendix B. 13

3.1.2

Exacte oplossing van Black-Scholes vergelijking

Voor de Black-Scholes vergelijking (9) bestaat er een unieke oplossing, die exact te berekenen is. We zullen de afleiding van deze oplossing niet presenteren, maar laten wel de resultaten zien. Voor een call optie C(S, t), met de eindvoorwaarde (10) en randvoorwaardes (11) en (12) geldt de oplossing C(S, t) = SN (d1 ) − Ke−r(T −t) N (d2 ),

(16)

waarbij N (x) de verdelingsfunctie is van een N (0, 1) stochast (de standaard normale verdeling) in x en

d1 = d2 =

ln(S/K) + (r + 12 σ 2 )(T − t) √ , σ T −t ln(S/K) + (r − 21 σ 2 )(T − t) √ . σ T −t

(17)

Voor een put optie P (S, t) met de eindvoorwaarde (13) en randvoorwaarden (14) en (15) geldt de oplossing P (S, t) = Ee−r(T −t) (1 − N (d2 )) + S(N (d1 ) − 1).

(18)

De exacte oplossing voor een call optie (16) is ge¨ımplementeerd in MATLAB, zie secties D.2 en D.3 van de MATLAB Appendix.

3.2 3.2.1

Heston model Afleiding van Heston optie vergelijking

In het Heston model beschouwen we de volatiliteit als een stochast, waardoor we model (2) krijgen voor een aandeel. Uit Itô’s lemma1 voor een functie V (S1 , S2 , t) van twee stochasten S1 en S2 en de tijd t volgt voor dV :

dV =

∂V ∂V ∂V 1 ∂2V ∂2V 1 ∂2V dt + dS1 + dS2 + B12 2 dt + B1 B2 dt + B22 2 dt ∂t ∂S1 ∂S2 2 ∂S1 ∂S2 2 ∂S1 ∂S2

(19)

waarbij B1 en B2 functies zijn als in (65). We beschouwen nu een porfolio Π bestaande uit een optie met waarde V (S, ν, t), −∆ van het onderliggende aandeel en −∆1 van een andere optie met waarde V1 (S, ν, t). De waarde van dit portfolio is nu 1

We gebruiken hier een uitbreiding naar twee stochasten op Itˆ o’s lemma uit A.1.

14

Π = V − ∆S − ∆1 V1 .

(20)

Nu kunnen we (19) invullen met S2 = ν en termen verzamelen zodat

1 2 ∂2V ∂2V 1 2 ∂2V ∂V dt + νS + ργSν + νγ dΠ = ∂t 2 ∂S 2 ∂S∂ν 2 ∂ν 2 ∂V1 1 2 ∂ 2 V1 ∂ 2 V1 1 2 ∂ 2 V1 −∆1 dt + νS + ργSν + νγ ∂t 2 ∂S 2 ∂S∂ν 2 ∂ν 2 ∂V ∂V1 ∂V ∂V1 dν. − ∆1 − ∆ dS + − ∆1 ∂S ∂S ∂ν ∂ν

(21)

We zien nu dat we ∆ en ∆1 zodanig kunnen kiezen, dat we een risicovrij portfolio krijgen, analoog aan bij de Black-Scholes analyse. De vergelijkingen waar ∆ en ∆1 aan moeten voldoen volgen uit de onderste regel van (21)

∂V ∂V1 − ∆1 −∆ = 0 ∂S ∂S ∂V ∂V1 − ∆1 = 0. ∂ν ∂ν Uit dezelfde redenering als bij de Black-Scholes analyse volgt dat dΠ = rΠdt, waardoor:

∂V 1 ∂2V ∂2V 1 ∂2V + νS 2 2 + ργSν + νγ 2 2 dt ∂t 2 ∂S ∂S∂ν 2 ∂ν 2 2 ∂V1 1 2 ∂ V1 ∂ V1 1 2 ∂ 2 V1 −∆1 + νS + ργSν + νγ dt ∂t 2 ∂S 2 ∂S∂ν 2 ∂ν 2 = rΠdt = r(V − ∆S − ∆1 V1 )dt.

dΠ =

(22)

We kunnen (22) nog herschrijven, zodat we vinden

∂V ∂t

=

∂V1 ∂t

2

2

2

∂ V + 12 νS 2 ∂∂SV2 + ργSν ∂S∂ν + 21 νγ 2 ∂∂νV2 + rS ∂V ∂S − rV 2

+ 12 νS 2 ∂∂SV21 +

∂V ∂ν 2 ∂ 2 V1 ργSν ∂S∂ν + 12 νγ 2 ∂∂νV21 ∂V1 ∂ν

1 + rS ∂V ∂S − rV

.

(23)

Hieruit volgt dat de linker- en rechterkant van deze vergelijking gelijk moeten zijn aan een functie g die alleen afhangt van S, ν en t. We kiezen g als g = κ(ν − ν) − θν. 15

(24)

Deze keuze wordt toegelicht in sectie 4.2. Gelijkstellen van (24) aan de linkerzijde van (23) en substitueren van κ0 = κ − θ en κ0 ν 0 = κν, levert ∂V 1 ∂2V ∂2V 1 ∂2V ∂V ∂V + νS 2 2 + ργSν + νγ 2 2 + rS − rV = κ0 (ν − ν 0 ) . ∂t 2 ∂S ∂S∂ν 2 ∂ν ∂S ∂ν

(25)

κ0 en ν 0 heten de risico-aangepaste parameters en worden gebruikt om de notatie te versimpelen. (25) is de optievergelijking voor het Heston model. Voor het oplossen van deze vergelijking zullen we randvoorwaarden nodig hebben. We zullen hier alleen de randvoorwaarden voor een call optie C(S, ν, t) beschouwen (de voorwaarden voor een put kunnen analoog worden gevonden, zie sectie 3.1.1). De eerste drie voorwaarden kunnen we direct halen uit de voorwaarden voor de Black Scholes vergelijking (zie ook [14]). De voorwaarde C(Smax , ν, t) = Smax voor Smax >> K wordt vaak herschreven als ∂C ∂S (Smax , ν, t) = 1.

C(S, ν, T ) = max(S(T ) − K, 0), C(0, ν, t) = 0, ∂C (Smax , ν, t) = 1. ∂S

(26)

Omdat we ook nog tweede orde afgeleides in ν hebben, zijn daar ook twee randvoorwaarden voor nodig. Over de keuze daarvoor lijkt geen uniek antwoord te vinden (zie bijvoorbeeld de discussie in [15]). In Heston’s oorspronkelijke paper worden de volgende twee randvoorwaarden voorgesteld.

∂C ∂C ∂C (S, 0, t) + κν (S, 0, t) − rC(S, 0, t) + (S, 0, t) = 0, ∂S ∂ν ∂t C(S, ∞, t) = S. rS

(27)

De voorwaarde op ν = 0 is simpelweg de resulterende vergelijking na invullen van ν = 0. Bij de numerieke experimenten blijkt dit echter geen praktische resultaten op te leveren. Daarom zullen wij gebruik maken van de voorwaarde C(S, 0, t) = max(S − K, 0).

(28)

De gedachte hierachter is dat als een aandeel geen volatiliteit heeft, de waarde van een optie constant blijft. Omdat we weten dat de waarde van de optie op t = T gelijk is aan max(S −K, 0), blijft deze waarde overal hieraan gelijk.

3.3

Heston-Hull-White model

In het Heston-Hull-White model beschouwen we de rente en de volatiliteit als stochasten, waardoor we het stelsel stochastische differentiaalvergelijkingen krijgen in (4). Voordat we aan de afleiding van een differentiaalvergelijking voor de optieprijs kunnen beginnen, zullen we eerst Itô’s lemma moeten uitbreiden. Dit wordt gedaan in de appendix. 16

3.3.1

Waarom afleiding van Heston-Hull-White optie vergelijking m.b.v. portfolio argumenten niet werkt

Bij de afleiding van de Black-Scholes en Heston optievergelijkingen hebben we telkens een portfolio argument gebruikt. In het geval van Heston is er nog een scheiding van variabelen toegepast om tot een vergelijking zonder stochastische termen te komen. In deze sectie zullen we laten zien waarom dit bij het Heston-Hull-White model niet met een vergelijkbaar portfolio kan. We gaan wederom een portfolio Π beschouwen, bestaande uit een optie met waarde V (S, ν, r, t), −∆ van het onderliggende aandeel, −∆1 van een andere optie met waarde V1 (S, ν, r, t) en −∆2 van een derde optie met waarde V2 (S, ν, r, t): Π = V − ∆S − ∆1 V1 − ∆2 V2 We maken gebruik van Itˆ o’s lemma voor 3 variabelen (70), waarbij we de volgende notatie introduceren om de afleiding leesbaarder te maken: Alleen de delen die dS, dν en dr bevatten zijn nog stochastisch, daarom schrijven we

dV

∂V ∂V ∂V ∂V dS + dν + dr + dt ∂S ∂ν ∂r ∂t 2 2 2 1 2∂ V 2∂ V 2∂ V + νS + νγ +η dt 2 ∂S 2 ∂ν 2 ∂r2 √ √ ∂2V ∂2V ∂2V + ρ13 νηS + ρ23 νγη dt + ρ12 νγS ∂S∂ν ∂S∂r ∂ν∂r ∂V ∂V ∂V dS + dν + dr = ∂S ∂ν ∂r + dV . =

Hierin bestaat dV uit alle niet-stochastische elementen van dV . Toepassen van de nieuwe notatie en ordenen van de termen bij een kleine verandering in het portfolio dΠ levert

dΠ = dV − ∆dS − ∆1 dV1 − ∆2 dV2 = dV − ∆1 dV 1 − ∆2 dV 2 ∂V ∂V1 ∂V2 + − ∆1 − ∆2 − ∆ dS ∂S ∂S ∂S ∂V1 ∂V2 ∂V − ∆1 − ∆2 dν + ∂ν ∂ν ∂ν ∂V ∂V1 ∂V2 + − ∆1 − ∆2 dr. ∂r ∂r ∂r

(29)

We zien nu dat we (29) risicovrij kunnen maken door ∆, ∆1 en ∆2 zo te kiezen dat de uitdrukking voor dS, dν en dr verdwijnen: 17

∂V ∂V1 ∂V2 − ∆1 − ∆2 −∆ = 0 ∂S ∂S ∂S ∂V1 ∂V2 ∂V − ∆1 − ∆2 = 0 ∂ν ∂ν ∂ν ∂V ∂V1 ∂V2 − ∆1 − ∆2 = 0. ∂r ∂r ∂r

(30)

Een oplossing voor het stelsel vergelijkingen in (30) is

∂V

∂V∂ν

∆2 = ∂V2 ∂ν

∆1 = ∆ =

∂V ∂r

2 − ∂V ∂r ∂r ∂V1 ∂r

−

∂V1 ∂ν

2 − ∆2 ∂V ∂r

∂V1 ∂r

∂V ∂V2 ∂V1 − ∆2 − ∆1 . ∂S ∂S ∂S

(31)

waarbij ∆1 en ∆ niet verder uitgewerkt zijn omdat de resultaten dan moeilijk leesbaar worden. Duidelijk is in ieder geval dat substitueren van (31) in (29) een niet-stochastische PDE geeft. Als we vervolgens gebruik maken van dΠ = rΠdt = r(V − ∆S − ∆1 V1 − ∆2 V2 )dt , vinden we r(V − ∆S − ∆1 V1 − ∆2 V2 )dt = dV − ∆1 dV 1 − ∆2 dV 2 .

(32)

Op deze vergelijking willen we tweemaal scheiding van variabelen toepassen (zoals bij het Heston model), om een exacte uitdrukking te vinden. Nu is het duidelijk waarom deze methode niet werkt bij het Heston-Hull-White model. Om scheiding van variabelen toe te passen, willen we (32) splitsen in een deel dat een functie is van V en een deel dat een functie is van V1 en V2 . Vanwege de gevonden uitdrukkingen voor ∆, ∆1 en ∆2 is dit niet mogelijk.

3.3.2

Heston-Hull-White optie vergelijking

We hebben gezien dat het portfolio argument bij het Heston-Hull-White model niet kan worden toegepast. Er zijn echter andere methodes ontwikkeld (zie [7]) die met behulp van Martingalen theorie de vergelijking hebben opgesteld. De afleiding hiervan ligt buiten het bereik van dit verslag, maar het uiteindelijke resultaat tonen we wel. De Heston-Hull-White vergelijking is

∂V ∂V ∂V ∂V + rS + λ(r − r) + κ(ν − ν) − rV ∂t ∂S ∂r ∂ν √ ∂2V √ ∂2V 1 ∂2V +ρ12 Sνγ + ρ13 Sη ν + ρ23 ηγ ν + νS 2 ∂S∂ν ∂S∂r 2 ∂r∂ν 2 2 2 ∂ V ∂ V 1 1 ∂ V 1 1 + νS 2 2 + η 2 2 + γ 2 ν + νS 2 2 = 0. 2 ∂S 2 ∂r 2 2 ∂ν 18

(33)

4 4.1

Details van het Heston model Motivatie voor Heston stochastische volatiliteit model

Het ligt voor de hand om de volatiliteit als een stochast te beschouwen (in tegenstelling tot een bekende reële functie zoals bij het traditionele Black-Scholes model), omdat praktische data aangeeft dat volatiliteit variabel en onvoorspelbaar is. Hieruit volgt echter niet direct waarom de vorm zoals gebruikt in het Heston model gekozen is. De verklaring hiervoor is de zogenaamde volatility smile die geobserveerd wordt in optiemarkten[10]. We zullen hier kort deze uitdrukking toelichten. We kunnen uit gegevens uit de markt een prijs Cmarkt vinden voor (bijvoorbeeld) een call √ met strike K en expiratietijd T . Dan kunnen we ons afvragen welke volatiliteit ν in de BlackScholes formule CBS (toegelicht in sectie 3.1.1 ) gesubstitueerd dient te worden om de marktprijs √ te krijgen. We willen dus de volgende vergelijking oplossen voor ν: √ CBS ( ν) = Cmarkt .

(34)

√ De resulterende ν noemen we de ge¨ımpliceerde volatiliteit. In [10] wordt aangetoond dat we deze vergelijking op kunnen lossen als max(S − Ke−r(T −t) , 0) ≤ Cmarkt ≤ S. √ De resulterende ν, blijkt onafhankelijk van K. M.a.w., het oorspronkelijke Black-Scholes model impliceert een volatiliteit die niet afhankelijk is van de hoogte van de exercise prijs. In de praktijk blijken deze twee echter wel van elkaar afhankelijk. Deze relatie wordt de ’volatilitysmile’ genoemd, naar aanleiding van de plot tussen volatiliteit en exercise prijzen (zie figuur 2). Het Heston model leidt tot vergelijkbare ’volatility smiles’, waardoor het dus beter bij de markt past.

4.2

Keuze van g(S, ν, t)

De keuze van g(S, ν, t) in (24) vereist nog enige toelichting. Wiskundig gezien zijn we vrij om elke functie g van S, ν en t te kiezen. Zonder algemeenheid te verliezen (zie [12]), kunnen we g kiezen als g(S, ν, t) = κ0 (ν − ν 0 ) − φ.

(35)

Hierin wordt φ(S, ν, t) de marktprijs van volatiliteitsrisico genoemd. Het is een functie die beschrijft hoe veel van de return van V be¨ınvloed wordt door het risico van ν. In [13] wordt vervolgens beschreven dat er economische argumenten zijn om te stellen dat deze φ proportioneel is aan ν. Daarom kiezen we φ = θν, 19

(36)

Figuur 2: De ’volatility smile’, zoals wordt geobserveerd in de markt. waardoor we uiteindelijk vinden: g = κ(ν − ν) − θν.

4.3

(37)

Coo ¨rdinatentransformatie

Bij het toepassen van numerieke methoden op partiële differentiaalvergelijkingen kan het zinvol zijn een co¨ ordinatentransformatie toe te passen op de variabelen. Een veelgebruikte transformatie bij het Heston model is om ln(S) te beschouwen in plaats van S (zie bijvoorbeeld [19]). In veel gevallen wordt dit gebruikt voor Monte Carlo simulatie. We kijken of de transformatie ook kan helpen bij onze numerieke schema’s. Het doel hierbij is een transformatie te vinden die de vergelijking vereenvoudigt, bijvoorbeeld als een aantal termen elkaar te laten cancellen. Voor de numerieke schema’s die wij zullen gebruiken is het vooral interessant te kijken wat er gebeurt met de randen van het gebied. Ter inleiding zullen we de transformatie S → ln(S), ν → ln(ν) bekijken. Daarna analyseren we een algemene coördinatentransformatie.

4.3.1

Co¨ ordinatentransformatie S → ln(S), ν → ln(ν)

We proberen de vergelijking in (25) te versimpelen door de volgende coördinatentransformatie toe te passen:

20

ˆ S → ln(S) = S, ν → ln(ν) = νˆ. Met behulp van de kettingregel kunnen we nu de partiële afgeleides uit (25) in S en ν omschrijven naar afgeleides in Sˆ en νˆ.

∂V ∂S

= =

∂2V ∂S 2

= = = = = = = = =

∂V ∂ν ∂2V ∂ν 2 ∂2V ∂ν∂S

= = =

∂V ∂ Sˆ ∂ Sˆ ∂S 1 ∂V S ∂ Sˆ ∂ 1 ∂V ∂S S ∂ Sˆ 1 ∂ ∂V 1 ∂V + − 2 ˆ S ∂S S ∂S ∂ Sˆ 1 ∂V 1 ∂ ∂V − 2 + ˆ S ∂S S ∂ Sˆ ∂S 1 ∂V 1 ∂ 1 ∂V − 2 + S ∂ Sˆ S ∂ Sˆ S ∂ Sˆ 1 ∂V 1 ∂V ∂ 1 1 ∂2V − 2 + + S ∂ Sˆ S ∂ Sˆ ∂ Sˆ S S ∂ Sˆ2 1 ∂V ∂ 1 ∂S 1 ∂2V 1 ∂V + + − 2 S ∂ Sˆ S ∂ Sˆ ∂S S ∂ Sˆ S ∂ Sˆ2 2 1 ∂V 1 1∂ V 1 ∂V + − 2S + − 2 S ∂ Sˆ S ∂ Sˆ S S ∂ Sˆ2 2 1 ∂V 1 1 ∂ V ∂V − 2 + − 2 ˆ ˆ S ∂S S S ∂S ∂ Sˆ 2 ∂V 1 ∂ V −2 S 2 ∂ Sˆ2 ∂ Sˆ 1 ∂V ν ∂ νˆ 1 ∂2V ∂V −2 ν 2 ∂ νˆ2 ∂ νˆ 2 1 ∂ V νS ∂ νˆ∂ Sˆ

Substitueren van bovenstaande transformaties in (25) levert de volgende vergelijking:

∂V 1 + ν ∂t 2

∂2V ∂V ν 0 ∂V ∂V ∂2V 1 γ2 ∂2V ∂V 0 +ργ + −2 +r −rV = κ 1 − . (38) −2 ˆ νˆ 2 ν ∂ νˆ2 ∂ νˆ ν ∂ νˆ ∂ Sˆ2 ∂ Sˆ ∂ S∂ ∂ Sˆ

ˆ Nu gebruiken we S = eS en ν = eνˆ , zodat we de volgende formule vinden voor Sˆ = ln(S) en νˆ = ln(ν):

21

∂V 1 νˆ + e ∂t 2

∂V ∂2V 1 2 −ˆν ∂ 2 V ∂V ∂V ∂2V 0 0 −ˆ ν ∂V −2 +ργ + γ e − 2 +r −rV = κ 1 − ν e . ˆ νˆ 2 ∂ νˆ2 ∂ νˆ ∂ νˆ ∂ Sˆ2 ∂ Sˆ ∂ S∂ ∂ Sˆ (39)

We kunnen ook de randvoorwaarden voor een call optie C(S, ν, t) uit (26) en (27) transformeren. We beginnen door op te merken dat S = 0 equivalent is met Sˆ = −∞ en S = ∞ met Sˆ = ∞. 1 ∂C Hetzelfde resultaat geldt uiteraard voor ν en νˆ. Gebruikmakend van ∂C ∂S = eSˆ ∂ Sˆ vinden we dan voor (26):

ˆ νˆ, T ) = max(S(T ˆ ) − K, 0), C(S, C(−∞, νˆ, t) = 0, ∂C ˆ (∞, νˆ, t) = eS . ∂ Sˆ

(40)

Voor (27) vinden we

ˆ

ˆ ∞, t) = eS , C(S, ∂C ˆ 1 ∂C ˆ ˆ −∞, t) + ∂C (S, ˆ −∞, t) = 0. r (S, −∞, t) + κν vˆ (S, −∞, t) − rC(S, ∂t e ∂ νˆ ∂ Sˆ

(41)

We zien dat in de laatste regel de uitdrukking e1νˆ voorkomt. Als νˆ → −∞, dan zal deze uitdrukking naar ∞ gaan, wat een probleem veroorzaakt.

4.3.2

Algemene co¨ ordinaten transformatie

Bij de transformatie S → ln(S), ν → ln(ν) hebben we gezien dat er geen numerieke voordelen ontstaan in de resulterende vergelijking en randvoorwaarden en dus is het dit geen goede transformatie. In deze sectie beschouwen we een algemene transformatie en beschrijven de eisen aan die transformatie om numerieke problemen te voorkomen. We laten eerst de gevolgen van algemene transformaties op de afgeleides van V zien. We gebruiken de volgende notatie

Sˆ = f (S), ˆ S = f −1 (S), dSˆ dS dS dSˆ d2 Sˆ dS 2

= f 0 (S), = [f −1 ]0 (S), = f 00 (S). 22

Nu zien we voor

∂V ∂S :

∂V ∂S

en voor

∂V ∂ Sˆ ∂ Sˆ ∂S ∂V = f 0 (S) ∂ Sˆ =

∂2V ∂S 2

∂2V ∂S 2

= = = = = = = =

∂ ∂V ∂S ∂S ∂V ∂ 0 f (S) ∂S ∂ Sˆ ∂V ∂ 0 ∂ ∂V 0 f (S) + f (S) ∂S ∂ Sˆ ∂ Sˆ ∂S ∂ ∂V ∂V + f 0 (S) f 00 (S) ˆ ∂S ∂S ∂ Sˆ ∂V ∂ ∂V 00 0 0 f (S) + f (S) f (S) ∂ Sˆ ∂ Sˆ ∂ Sˆ ∂V ∂V ∂ 0 ∂2V 00 0 0 f (S) + f (S) f (S) + f (S) ∂ Sˆ ∂ Sˆ ∂ Sˆ ∂ Sˆ2 ∂V 00 ∂2V ∂V 0 −1 0 ˆ 0 00 + f (S) f (S)[f ] (S) + f (S) f (S) ∂ Sˆ ∂ Sˆ ∂ Sˆ2 2 ∂ V ˆ ∂V . (f 0 (S))2 + f 00 (S) 1 + f 0 (S)[f −1 ]0 (S) ∂ Sˆ2 ∂ Sˆ

ˆ toe te passen. We kunnen de transformaties als functies van Sˆ beschouwen door S = f −1 (S) Dezelfde analyse kunnen we toepassen voor transformaties in ν. Nu kunnen we (42) toepassen op de oorspronkelijke Heston vergelijking, om de getransformeerde formule te vinden. We passen Sˆ = f (S) en νˆ = u(ν) toe. We vinden dan de volgende vergelijking

∂V 1 ∂2V ˆ ∂V + νS 2 (f 0 (S))2 + f 00 (S) 1 + f 0 (S)[f −1 ]0 (S) ∂t 2 ∂ Sˆ2 ∂ Sˆ 2 2 ∂ V 1 ∂ V + ργSνf 0 (S)u0 (ν) + νγ 2 (u0 (v))2 2 ˆ νˆ 2 ∂ νˆ ∂ S∂ ∂V ∂V + u00 (ν) 1 + u0 (ν)[u−1 ]0 (ˆ ν) + rS − rV ∂ νˆ ∂ Sˆ ∂V = κ0 (ν − ν 0 )u0 (ν) . ∂ νˆ

(42)

We gaan nu bekijken hoe we de problemen die we bij de transformatie S → ln(S), ν → ln(ν) tegenkwamen kunnen generaliseren. Er zijn twee type problemen:

23

1. De uitdrukkingen voor de afgeleides worden snel te groot. De reden dat dit een probleem vormt, is dat een specifieke afgeleide gaat overheersen in de uiteindelijke oplossing. We zoeken dus een transformatie, waardoor de termen voor de afgeleides begrensd zijn op de randen van het getransformeerde gebied. 2. De transformatie zorgt voor een extra truncatie in het numerieke schema. In de oorspronkelijke co¨ ordinaten S en ν kiezen bekijken we het gebied [0, ∞]×[0, ∞]. Voor het numerieke schema zullen we dit moeten trunceren tot [0, Smax ] × [0, νmax ]. Als door de transformatie de rand S = 0 of ν = 0 op −∞ komt, moeten we een extra truncatie toevoegen. Hierdoor wordt de oplossing minder nauwkeurig. Als we de S termen in (43) bekijken, zien we dat er twee uitdrukkingen voor de afgeleides staan die problemen op kunnen leveren. Om probleem 1 te voorkomen, moeten we f zo kiezen dat Sf 0 (S) en S 2 (f 0 (S))2 niet te snel groot worden. Hieruit volgt dat f 0 (S) = S1 (oftewel f (S) = ln(S)) een logische kandidaat is, hierdoor maken we de oorspronkelijke vergelijking ook nog gemakkelijker omdat de S en S 2 term vervallen. 0 −1 0 00 ˆ Er ontstaat echter dan een probleem in de uitdrukking f (S) 1 + f (S)[f ] (S) . Hier introduceren we dan namelijk de uitdrukking f 00 (S) = − S12 , die niet is gedefinieerd op de rand S = 0. Voor de transformatie in ν is er zelfs geen geschikte kandidaat te vinden. Deze moet dan de uitdrukking νu0 (ν) en ν(u0 (ν))2 versimpelen, zonder een singulariteit te creëren op de rand ν = 0. We zien dat zo’n transformatie niet bestaat. Daarom kunnen we concluderen dat het transformeren van het stelsel niet zinvol is.

4.4

Feller conditie

Om te garanderen dat de variantie in (2) positief blijft (ν > 0), moet aan de zogenaamde Feller conditie worden voldaan [5]: 2κν ≥ γ 2

(43)

In [5] wordt geschreven dat het proces 2ζνt voor 0 < s < t beschreven wordt door een non-centrale chi-kwadraat verdeling met graad q en non-centraliteit parameter 2ζνs e−κ(t−s) , χ2 (q, 2ζνs e−κ(t−s) ) waarbij:

ζ = q =

2κ e−κ(t−s) )γ 2

(1 − 2κν −1 γ2

De Feller conditie is dus equivalent met q ≥ 0. In de praktijk blijkt vaak niet aan de Feller conditie te worden voldaan. We zullen bij de numerieke resultaten zien dat de fouten inderdaad groter worden wanneer niet voldaan is aan de Feller conditie. 24

5 5.1

Numerieke implementatie Eindige differenties voor Black-Scholes

Voor het oplossen van de Black Scholes PDV (9) gaan we eindige differenties gebruiken. We volgen de aanpak uit hoofdstuk 24 van [1], waarbij we voorwaartse differentie in de tijd t gebruiken en centrale differentie in de waarde van het onderliggende aandeel S. Hiervoor passen we de co¨ ordinatentransformatie τ = T − t toe, zodat we onze eindtijdconditie om kunnen schrijven naar een begintijdconditie. Na deze transformatie vinden we ∂V 1 ∂2V ∂V − σ 2 S 2 2 − rS + rV = 0. ∂τ 2 ∂S ∂S

(44)

We gaan de optiewaarde V bepalen op een NS × Nτ grid in S en τ . We discretizeren S en τ als

Sj

= jh,

τi = ik. Hierin zijn h en k de stapgroottes in respectievelijk S en t. Gebruikmakend van de voorwaartse en centrale differenties, vinden we voor de partiële afgeleides op het punt (jh, ik)

∂V ∂τ ∂V ∂S ∂2V ∂S 2

= = =

Vji+1 − Vji

, k i i Vj+1 − Vj−1 , 2h i i Vj+1 − 2Vji + Vj−1 . h2

Nu vinden we voor (44) Vji+1 − Vji k

i i i − Vj−1 V i − 2Vji + Vj−1 Vj+1 1 2 2 j+1 − σ (jh) − rjh + rVji = 0. 2 h2 2h

(45)

Dit kunnen we omschrijven naar een uitdrukking voor Vji+1 1 1 i i i i Vji+1 = Vji (1 − rk) + kσ 2 j 2 (Vj+1 − 2Vji + Vj−1 ) + rkj(Vj+1 − Vj−1 ). 2 2 Nu kunnen we een matrix-vector representatie voor deze vergelijking opstellen. Laat

25

(46)

V1i V2i .. .

   Vi =  

VNi S −1

    ∈ RNS −1 

(47)

de numerieke oplossing zijn op τi . Voordat we verder gaan, moeten we eerst kijken naar welke randvoorwaarden we hebben. Als we een call optie beschouwen, hebben we de randvoorwaarden V (0, τ ) = 0 en C(Smax , τ ) = Smax . Dit kunnen we omschrijven naar

V0i = 0 i VN S

= Smax

(48)

de numerieke oplossing zijn op τi . Verder hebben we de beginconditie (τ = 0 geeft t = T ) V (S, τ = 0) = max(S(τ = 0) − K, 0), die we om kunnen schrijven naar Vj0 = max(Sj − K, 0) = max(jh − K, 0)

(49)

We gaan nu op zoek naar een Matrix-Vector representatie van de vorm Vi+1 = F Vi + pi ,

(50)

waarbij pi gebruikt wordt om de randvoorwaarden te verwerken. Uit (46) volgt dat voor F geldt 1 1 F = (1 − rk)I + kσ 2 D2 T2 + krD1 T1 2 2 waarbij

26

(51)



D1

D2

T1

T2

 ... ... 0   .. ..  0 2  . 0 .    ..  .. .. =  . 0 , . 3 .    .. . . . . . .   .  . . . 0 0 . . . . . . 0 NS − 1  2  1 0 ... ... 0   ..  0 22 0 . . .  .     .. .. 2 =  ... , . 0 3 .    .. . . . . . .   .  . . . 0 2 0 . . . . . . 0 (NS − 1)   0 1 0 ... ... 0  .  .. ..  −1 0 . . ..  1    ..  . . . . . . . .  0 . .  . . . , =    .. . . . . . . . .  . . 0  . . .     .. . . . .  . . −1 0 1  . 0 . . . . . . 0 −1 0   −2 1 0 ... ... 0  ..   1 −2 1 . . . . . . .     ..  .. .. ..   0 . . . . 1 . =    .. . . . . . . . .  . . 0  . . .     .. .. ..  . . 1 −2 1  . 0 ... ... 0 1 −2 1

0

De keuze van deze matrices volgt direct uit (46) op de volgende manier 1 1 i i i i − 2Vji + Vj−1 ) + rk j (Vj+1 − Vj−1 ), Vji+1 = Vji (1 − rk) + kσ 2 j 2 (Vj+1 |{z} |{z} 2 2 |{z} | {z } | {z } D2

I

D1

T2

T1

waarbij de underbraces aangeven in de rij van welke matrix dat gedeelte van de formule terugkomt. We zien dat T1 en T2 op de eerste een laatste rij niet overeenkomen met (46). Dit kunnen we oplossen met de vector pi . Voor j = 1 vinden we

1 1 V1i+1 = V1i (1 − rk) + kσ 2 · 12 · (V2i − 2V1i + V0i ) + rk · 1 · (V2i − V0i ) 2 2 1 2 2 1 1 1 i i i 1 (V2 − 2V1 ) + rk |{z} 1 V2i + kσ 2 V0i − rkV0i , = V1 (1 − rk) + kσ |{z} |{z} |{z} |2 | {z } 2 2 2 {z } D2 D1 I

T2

27

T1

pi

waarbij de underbraces dezelfde betekenis hebben als voorheen. Voor j = NS − 1 vinden we

1 1 VNi+1 = VNi S −1 (1 − rk) + kσ 2 (NS − 1)2 (VNi S − 2VNi S −1 + VNi S −2 ) + rk(NS − 1)(VNi S − VNi S −2 ) S −1 2 2 1 1 = VNi S −1 (1 − rk) + kσ 2 (NS − 1)2 (−2VNi S −1 + VNi S −2 ) + rk (NS − 1) −VNi S −2 | {z } | 2 {z } 2 | {z } | {z } | {z } D2

I

D1

T2

T1

1 1 + kσ 2 (NS − 1)2 VNi S + rk(NS − 1)VNi S . {z 2 } |2 pi

Hieruit volgt voor pi      pi =     

1 2 2 k(σ

− r)V0i 0 .. . .. . 0 1 2 i 2 k(NS − 1)(σ (NS − 1) + r)VNS

     .    

Het schema Vi+1 = F Vi + pi kan nu worden ge¨ımplementeerd in MATLAB. Dit is voor een call-optie gedaan, zie sectie D.1 van de MATLAB appendix.

5.2

Eindige differenties voor Heston

We gaan eindige differenties toepassen op de ongetransformeerde Heston vergelijking 1 ∂2V ∂2V 1 ∂2V ∂V ∂V ∂V + νS 2 2 + ργSν + νγ 2 2 + rS − rV = κ0 (ν − ν 0 ) . ∂t 2 ∂S ∂S∂ν 2 ∂ν ∂S ∂ν

(52)

We maken gebruik van voorwaartse differentie in de tijd t en centrale differenties in S en ν. Voor de convectie termen zullen we upwind discretisatie gebruiken, zie appendix C. We passen eerst de transformatie τ = T − t toe, zodat we de eindtijdconditie om kunnen schrijven naar een beginwaarde. We schrijven direct om naar een uitdrukking voor ∂V ∂τ , zodat we vinden ∂V 1 ∂2V ∂2V 1 ∂2V ∂V ∂V = νS 2 2 + ργSν + νγ 2 2 + rS − rV + κ0 (ν − ν 0 ) . ∂τ 2 ∂S ∂S∂ν 2 ∂ν ∂S ∂ν We passen de volgende discretisaties toe op τ, S en ν:

τi = ik, Sj

= jh,

νm = ml. 28

(53)

i voor V op Hierin zijn k, h en l de stapgroottes in τ, S en ν respectievelijk. We schrijven nu Vj,m τi , Sj en νm . Voor de partiële afgeleides in (τi , Sj , νm ) gebruiken we

∂V ∂τ ∂V ∂S ∂V ∂ν ∂2V ∂S 2 ∂2V ∂ν 2 ∂2V ∂S∂ν

= = = = = =

i+1 i Vj,m − Vj,m

k

,

(54)

i i Vj+1,m − Vj−1,m , 2h i i Vj,m+1 − Vj,m−1 , 2l i i +Vi Vj+1,m − 2Vj,m j−1,m , 2 h i i +Vi Vj,m+1 − 2Vj,m j,m+1 , l2 i i i i Vj+1,m+1 − Vj−1,m+1 − Vj+1,m−1 + Vj−1,m−1 . 4lh

(55) (56) (57) (58) (59)

Deze discretisaties passen we toe op (53). We versimpelen de notatie door de indices i, j en i m alleen te vermelden als ze afwijkend zijn (we schrijven bijvoorbeeld Vj+1 voor Vj+1,m ). Nu i+1 vinden we voor V

kmlj 2 ργkmj (Vj+1 − 2V + Vj−1 ) + (Vj+1,m+1 − Vj−1,m+1 − Vj+1,m−1 + Vj−1,m−1 ) 2 4 rkj γ 2 km (Vm+1 − 2V + Vm−1 ) + (Vj+1 − Vj−1 ) + 2l 2 kmκ0 kν 0 κ0 − (rk − 1)V − (Vm+1 − Vm−1 ) + (Vm+1 − Vm−1 ) 2 2l

V i+1 =

We willen dit omschrijven naar een matrix-vector notatie van de vorm Vi+1 = F Vi . We beginnen door ons [0, Nν ] × [0, NS ] grid horizontaal te nummeren in de ν-richting. Hierdoor wordt V een vector van lengte (Nν + 1)(NS + 1). Door de nummering in de ν-richting te kiezen, gebruiken we (ν = m, S = j) → m + 1 + j(Nν + 1) en vinden we

Vj,m = Vn , Vj,m+1 = Vn+1 , Vj+1,m = Vn+Nν +1 , Vj+1,m+1 = Vn+Nν +2 . waarbij n = m + 1 + j(Nν + 1). Nu kunnen we F schrijven als 29

F =

kl ργk γ2k rk kκ0 kν 0 κ0 D1 T1 + D2 T2 + D3 T3 + D4 T4 − (rk − 1)I − D5 T5 + T6 , 2 4 2l 2 2 2l

met de keuzes van D en T analoog aan die uit de behandeling van Black-Scholes, gebruikmakend van de nummering in de ν-richting. Implementatie van de randen is op dit moment eenvoudig, omdat we daar een Dirichlet voorwaarde hebben. We kunnen in de rijen van de elementen op die randen de corresponderende rij uit de identiteitsmatrix toepassen. Bij het behandelen van de resultaten zullen we echter zien dat we ook alternatieve randvoorwaarden gebruiken. Deze worden ge¨ımplementeerd met behulp van een upwind schema. In de appendix C wordt dit schema verder toegelicht.

6

Resultaten

In deze sectie analyseren we de resultaten van de numerieke schema’s ontwikkeld voor de BlackScholes en Heston optievergelijking. We vergelijken deze resultaten met de exacte oplossingen.

6.1

Black-Scholes model

Tenzij anders vermeld, wordt in deze sectie een call optie behandeld en zijn de volgende parameters gebruikt:

K = 100, σ = 0.2, r = 0.05, T

= 1,

Smax = 200. In figuur 3, zien we de resultaten voor NS = 100, Nt = 500. We zien direct dat er een fout is op de rand S = Smax . Deze fout wordt veroorzaakt door de truncatie van het domein. De keuze voor de randvoorwaarde V (Smax , t) = Smax is gebaseerd op het argument dat K verwaarloosbaar klein is t.o.v. Smax (Smax >> K). In het numerieke schema geeft dit dus een fout van grootte Smax − K op de rand. Deze fout werkt door in de rest van het schema, zoals te zien in figuur 4. Naast de fout veroorzaakt door de truncatie, wordt er ook een numerieke fout ge¨ıntroduceerd door het differentie-schema. De voorwaartse differentie in de tijd veroorzaakt een fout van O(k), de centrale differentie in S veroorzaakt een fout van O(h2 ). Uit figuur 4 volgt dat de fout veroorzaakt door de truncatie van het domein groter is dan de fout in de differenties. We kunnen op twee manieren met de domein-fout omgaan. De eerste mogelijkheid is om Smax groter te kiezen. Het probleem hierbij is dat we NS dan evenredig groter moeten maken, omdat we anders een grotere fout in de differentie introduceren. Dit vertraagt het schema.

30

Figuur 3: Plots van de oplossing van de Black-Scholes vergelijking. We zien dat de truncatie van gebied een fout in de numerieke oplossing veroorzaakt. Boven: Exacte oplossing. Onder: Numerieke oplossing. Een alternatieve oplossing vinden we door de randvoorwaarde anders te kiezen. In plaats van te gebruiken dat V (Smax , t) = Smax , kiezen we nu dat V (Smax , t) = Smax − K. Het theoretische argument dat limS→∞ V (S, t) = S blijft dan geldig. Praktisch voorkomt het de domein-fout die verooraakt wordt door een eindige Smax . Deze randvoorwaarde is toegepast, het resultaat (bij NS = 100, Nt = 500) is zichtbaar in figuur 5. De fout op t = 0 wordt afgebeeld in figuur 6. Tot slot moeten we nog een korte opmerking maken over de stabiliteit van het schema. Omdat we gebruik maken van een voorwaartse differentie in de tijd en een centrale differentie in S, geldt er een stabiliteitseis van de vorm

31

Figuur 4: Grafiek van de fout van het numerieke schema op t = 0 (τ = T ). We zien dat er een fout rond S = Smax is van orde S − K. Deze fout wordt veroorzaakt door de truncatie van het S-domein.

k < Ch2 .

(60)

We kunnen C theoretisch bepalen door de versterkingsfactor van ons schema te onderzoeken. Een alternatieve methode is te bekijken bij welke verhouding de numerieke tests onstabiele resultaten geven. Hieruit is gebleken dat we stabiliteit hebben voor C < 1/1500.

32

Figuur 5: Plots van de oplossing van de Black-Scholes vergelijking met de nieuwe randvoorwaarde op S = Smax . We zien dat er een veel betere overlap is tussen de waardes. Boven: Exacte oplossing. Onder: Numerieke oplossing.

33

Figuur 6: Grafiek van de fout van het numerieke schema op t = 0 (τ = T ) bij de nieuwe randvoorwaarde op S = Smax . We zien dat de fout veel kleiner is dan bij de oorspronkelijke randvoorwaarde. De resterende fout komt door de differentie.

34

6.2

Heston model

Het is belangrijk op te merken dat de plots nu allemaal zijn op het vaste tijdstip t = 0 (τ = T ). Het domein van de plots is dus een (S, ν) domein in tegenstelling tot het (S, t) domein bij de Black-Scholes analyse. Tenzij anders vermeld, worden in deze sectie call-opties behandeld en zijn de volgende waardes voor de parameters gebruikt:

K = 100, r = 0, ρ = −0.9, ν 0 = 0.04, γ = 1, κ0 = 5, NS = 200, Nv = 20, Nt = 10000, Smax = 200, νmax = 0.2, T

= 1.

We beginnen de analyse door een plot te bekijken met de randvoorwaardes zoals gegeven in sectie 3.2.1. In figuur 7 zijn de numerieke en exacte oplossing afgebeeld met de aangepaste randvoorwaarde op ν = 0. We zien direct dat we op de rand S = Smax dezelfde fout introduceren als bij het Black-Scholes model. De truncatie van het domein introduceert een fout van grootte S − K, omdat de achterliggende gedachte bij de voorwaarde V (Smax , ν, t) = Smax is dat K verwaarloosbaar klein is. Net als bij de Black-Scholes analyse kunnen we dit probleem oplossen door de voorwaarde te vervangen door V (Smax , ν, t) = Smax − K. We vinden dan een resultaat van de vorm als in grafiek 8. Vanaf nu gebruiken we dus de randvoorwaarde V (Smax , ν, t) = Smax − K.

(61)

Na het toepassen van de nieuwe randvoorwaarde op S = Smax , vinden we de plot in figuur 8. We zien dat er nog steeds een fout van orde S − K gemaakt wordt, maar nu alleen op de rand ν = νmax ! De randvoorwaarde daar is V (S, νmax , t) = S, maar in de literatuur is deze voorwaarde niet volledig geaccepteerd (zie sectie 3.2.1). Net als op de rand S = Smax , zullen we de voorwaarde vervangen door S − K. De randvoorwaarde op ν = νmax is nu dus V (S, νmax , t) = S − K.

35

(62)

Figuur 7: Plots van de oplossing van de Heston vergelijking. We zien dat de truncatie van gebied een fout in de numerieke oplossing veroorzaakt. Boven: Exacte oplossing. Onder: Numerieke oplossing.

36

Figuur 8: Numerieke oplossing met alternatieve randvoorwaarde op S = Smax . Het truncatieprobleem is op die rand verholpen, maar bestaat nog steeds op ν = νmax .

37

Na toepassing van de verbeterde randvoorwaardes, vinden we uiteindelijk een grafiek van de vorm die we zouden verwachten, zie figuur 9. Als we dit resultaat vergelijken met de exacte oplossing, kunnen we een plot van de fout maken. Dit is gedaan in figuur 10. Hierin is de absolute waarde van het verschil tussen de numerieke en de exacte oplossing afgebeeld. We zien dat er een opvallend piek loopt over het gebied rond S − K.

Figuur 9: Numerieke oplossing met alternatieve randvoorwaardes op S = Smax en ν = νmax . Het truncatieprobleem is nu verholpen. De piek blijkt inderdaad afhankelijk te zijn van de keuze van K. In figuur 11 is de fout afgebeeld voor K = 150, daarin zien we dat de piek is meegeschoven. De vraag is hoe we deze piek kunnen verklaren. Een logische gedachte is dat deze fout veroorzaakt wordt door de fouten bij het toepassen van differenties. Bij experimenten met zeer kleine stapgroottes bleef de fout echter bestaan en van dezelfde grootte2 . Dit is te zien in figuur 12, die gemaakt is met stapgroottes h = .25, l = 0.025, k = 1.6 · 10−5.

2

Bij deze stapgroottes is niet meer gebruik gemaakt van de Matrix-Vector vorm van de numerieke schema’s zoals beschreven in sectie 5. De benodigde matrices werden dan te groot voor MATLAB. In plaats daarvan werd een recursief schema gebruikt, zie de MATLAB appendix.

38

Figuur 10: Plot van de absolute waarde van het verschil tussen de numerieke en de exacte oplossing. We zien een opvallende piek in de fout rond het gebied S − K.

Figuur 11: Plot van de absolute waarde van het verschil tussen de numerieke en de exacte oplossing bij K = 150. We zien dat de piek in de fout rond het gebied S − K is meegeschoven.

39

Figuur 12: Plot van de absolute waarde van het verschil tussen de numerieke en de exacte oplossing bij h = .25, l = 0.025, k = 1.6e − 5. We zien dat de piek in de fout rond het gebied S − K nog steeds bestaat en even groot is gebleven.

40

We kunnen wel het foutenprofiel in het centrum afvlakken door νmax groter te kiezen. In figuur 13 zien we dit, hierbij is νmax = 1. Helaas blijft het probleem aan de rand bestaan. De oplossing hiervoor vinden we in [26], waar een alternatieve randvoorwaarde wordt opgesteld. Daar maken ze gebruik van de voorwaarde ∂V (S, νmax , t) = 0. ∂ν

(63)

Met deze voorwaarde is het probleem in de rand ν = νmax verholpen! In figuur 14 is de voorwaarde ge¨ımplementeerd met eenzijdige differenties, wederom met νmax = 1.

Figuur 13: Plot van de absolute waarde van het verschil tussen de numerieke en de exacte oplossing met νmax = 1. Hoewel het foutprofiel in het centrum afgevlakt is, blijft het probleem aan de rand bestaan.

41

Figuur 14: Plot van de absolute waarde van het verschil tussen de numerieke en de exacte oplossing met randvoorwaarde ∂V ∂ν (S, νmax , t) = 0. Het probleem aan de rand van ν = νmax is hiermee opgelost!

42

S=0 V =0

S = Smax V =S V = max(S − K, 0) ∂V ∂S = 0 ∂V ∂S = 1

ν=0 V =0

rS ∂V ∂S

+

∂V ∂ν ∂V κν ∂ν

=0 − rV +

∂V ∂t

ν = νmax V =S V = max(S − K, 0) ∂V ∂ν = 0 =0

Tabel 1: Overzicht van de verschillende gebruikte randvoorwaardes. De combinatie van de rode voorwaarden leveren het beste foutprofiel. Helaas is dit de laatste verbetering die we met behulp van randvoorwaarden toe hebben kunnen passen. Er zijn nog veel verschillende randvoorwaardes geprobeerd, in tabel 1 is een overzicht van deze voorwaarden te zien. De gebruikte voor de beste resultaten zijn hierin dikgedrukt. Het probleem op de rand S = Smax dat nog overblijft, is een gevolg van de code uit [14] die gebruikt wordt voor de exacte oplossing. Deze blijkt slecht te functioneren op de rand S = Smax voor grote T . Dit is te zien in figuur 15, waarin we het foutprofiel zien voor T = 0.01. We zien dat de fout op S = Smax verdwenen is. Merk op dat de schaal van het foutprofiel ook verkleind is, omdat we minder stappen bekijken. Alleen de fouten op ν = 0 en S = K blijven dus bestaan. De fout op S = K wordt vermoedelijk veroorzaakt door de ’kink’ in de beginvoorwaarde V = max(S − K, 0). De fout op ν = 0 is niet van groot belang, omdat deze situatie in praktijk nauwelijks voorkomt. De numerieke code is ook ge¨ımplementeerd voor put-opties, zie MATLAB appendix D. In de numerieke oplossing maken we dan gebruik van andere voorwaarden in S en de andere beginvoorwaarde(zoals gepresenteerd in sectie 3.2.1). Voor de exacte oplossing maken we gebruik van de put-call pariteit.

Figuur 15: Plot van de absolute waarde van het verschil tussen de numerieke en de exacte oplossing met T = 0.01. Het probleem aan de rand S = Smax is hiermee verholpen. 43

7

Conclusie

In dit verslag hebben we verschillende modellen voor de waardering van opties bekeken. We hebben voor Black-Scholes en Heston laten zien hoe met behulp van portfolio argumenten de stochastische differentiaal vergelijkingen omgeschreven kunnen worden tot partiële differentiaal vergelijkingen. Voor Heston-Hull-White is aangetoond waarom dit argument daar niet werkt. Vervolgens zijn op Black-Sholes en Heston eindige differenties toegepast. Hierbij is aangetoond dat en co¨ ordinatentransformatie geen numerieke voordelen oplevert. Bij de numerieke resultaten hebben we kunnen zien dat de randvoorwaardes voor de modellen niet goed praktisch ge¨ımplementeerd kunnen worden. Vanwege de truncatie van het domein ontstaan er fouten in de numerieke oplossingen. Er zijn alternatieve randvoorwaarden gepresenteerd die beter presteren en met dezelfde theoretische argumenten onderbouwd worden. Er blijft echter nog steeds een opvallende fout bestaan in het gebied rond S = K en op de rand ν = 0, die verder onderzocht kunnen worden.

8

Referenties

Referenties [1] D.J. Higham: An Introduction to Financial Option Valuation, Cambridge University Press (2004), ISBN 978-0-521-83884-9. [2] L.A. Grzelak & C.W. Oosterlee: On The Heston Model with Stochastich Interest Rates, ISSN 1389-6520. [3] http://en.wikipedia.org/wiki/Heston model [4] C. Kahl & P. J¨ ackel: Not-so-complex logarithms in the Heston model, Wilmott Magazine: 74-103. [5] F. Fang: The COS Method: An Efficient Fourier Method for Pricing Financial Derivatives, proefschrift TU Delft (2010). [6] L.A. Grzelak, C.W. Oosterlee, S. van Weeren: Extension of Stochastic Volatility Models with Hull-White Interest Rate Process, ISSN 1389-6520. [7] F.H.C. Naber: Fast solver for the three-factor Heston-Hull/White problem, MSc Thesis ING Amsterdam & TU Delft. [8] http://en.wikipedia.org/wiki/Hull%E2%80%93White model [9] P. Wilmot, J. Dewynne, S. Howison: Option Pricing: Mathematical Models and Computation, Oxford Financial Pss, ISBN 0 9522082 0 2 [10] M. Ng: MSc Thesis TU Delft, 2005 [11] http://en.wikipedia.org/wiki/Interest rate cap and floor

44

[12] J. Gatheral & M. Lynch: Stochastic Volatility and Local Volatility, Case Studies in Financial Modelling Course Notes, Courant Institute of Mathematical Sciences, Fall Term, 2004. [13] J.B. Wiggins: Option values under stochastic volatility, Journal of Financial Economics 19 (1987), 351-372. [14] S. Lin: Finite Difference Schemes for Heston Model, 06/2008.

[15] http://www.wilmott.com/messageview.cfm?catid=34&threadid=58758&FTVAR MSGDBTABLE=&ST [16] S.L. Heston: A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. [17] J. van Kan, A. Segal, F. Vermolen: Numerical Methods in Scientific Computing, VSSD, ISBN-13 978-90-71301-50-6. [18] C. Vuik, P. van Beek, F. Vermolen, J. van Kan: Numerieke Methoden vor Differentiaalvergelijkingen, VSSD, ISBN 978-90-65621-79-5. [19] L. Andersen: Efficient Simulation of the Heston Stochastic Volatility Model, Banc of America Securities. [20] R. Lord, R. Koekkoek, D. van Dijk: A comparison of biased simulation schemes for stochastic volatility models, Quantitative Finance (March 2008). [21] J.C. Cox, J.E. Ingersoll, S.A. Ross: A Theory of the Term Structure of Interest Rates, Econometrica: Volume 53, Issue 2 (March 1985), p. 385-408. [22] L.A. Grzelak, C.W. Oosterlee, S. van Weeren: Efficient Option Pricing with Multi-Factor Equity-Interest Rate Hybrid Models, ISSN 1389-6520. [23] W. Feller: Two Singular Diffusion Problems, The Annals of Mathematics, Second Series, Vol. 54, No. 1 (July 1951), p. 173-182. [24] C. Munk: Introduction to the Numerical Solution of Partial Differential Equations in Finance, October 2007. [25] D.B. Ntwiga: Numerical Methods for the Valuation of Financial Derivatives, MSc Thesis, November 2005. [26] C. O’ Sullivan, S. O’ Sullivan: Pricing Options under Heston’s Stochastic Volatility Model via Accelerated Explicit Finite Differencing Methods, June 2010. [27] G. Winkler, T. Apel, U. Wystup: Valuation of Options in Heston’s Stochastic Volatility Model Using Finite Element Methods, October 24 2001. [28] J. Wang: Convexity of option prices in the Heston model, January 2007.

A

Itˆ o’s lemma

Bij de afleiding van een Partiële Differentiaal Vergelijking voor de prijs van een optie, maken we vaak gebruik van Itˆ o’s lemma. Het lemma voor 1 of 2 variabele(n) komt in veel literatuur voor en wordt hieronder alleen geciteerd. Itô’s lemma voor 3 variabelen is niet uit externe bronnen gehaald en wordt daarom specifiek vermeld. 45

A.1

Itˆ o’s lemma voor 1 variabele

Itô’s lemma relateert een kleine verandering in een functie van een stochast met een kleine verandering in de stochast zelf. Voor een afleiding van de stelling, zie [9]: Voor een stochast G, die beschreven kan worden als dG = A(G, t)dW + B(G, t)dt

(64)

1 2 d2 f df df dt. df = A dW + B + A dG dG 2 dG2

(65)

geldt voor elke functie f (G):

Specifiek geldt voor elke functie f (S), met S als in (1):

df

df 1 d2 f (σSdW + rSdt) + σ 2 S 2 2 dt dS 2 dS df df 1 2 2 d2 f = σS dW + rS + σ S dt. dS dS 2 dS 2 =

(66)

We kunnen (66) uitbreiden voor een functie f (S, t) die ook afhankelijk is van de tijd, dan volgt uit Taylor reeks expansie

df = σS

A.2

∂f 1 ∂2f ∂f ∂f dW + rS + σ2S 2 2 + dt. ∂S ∂S 2 ∂S ∂t

(67)

Itˆ o’s lemma voor 3 variabelen

In deze sectie ontwikkelen we een uitbreiding van Itø’s lemma voor een kleine verandering in een functie V (S1 , S2 , S3 , t) van 3 stochasten en de tijd. Dit resultaat is niet uit externe bronnen gehaald en wordt daarom hier specifiek vermeld. De stochasten Si zijn van de vorm Si = ai dt + bi dWi

(68)

waarbij de onderlinge Wiener processen gecorreleerd zijn volgens dWi dWj = ρij dt. Door V te ontwikkelen als Taylor-polynoom in S1 , S2 , S3 en t en gebruik te maken van het feit dat dWi2 → dt als dt → 0, vinden we voor kleine veranderingen in dV :

dV

=

∂V ∂V ∂V ∂V dS1 + dS2 + dS3 + dt ∂S1 ∂S2 ∂S3 ∂t 2 1 ∂ V 2 ∂2V 2 ∂2V 2 + b + b + b dt 2 ∂S12 1 ∂S22 2 ∂S32 3 ∂2V ∂2V ∂2V + b1 b2 ρ12 + b1 b3 ρ13 + b2 b3 ρ23 dt ∂S1 ∂S2 ∂S1 ∂S3 ∂S2 ∂S3 46

(69)

Toegepast op het Heston-Hull-White model (4), levert Itô’s lemma (69) voor een optie (of andere functie) V (S, ν, r, t):

dV

B

=

∂V ∂V ∂V ∂V dS + dν + dr + dt ∂S ∂ν ∂r ∂t 2 2 2 1 2∂ V 2∂ V 2∂ V νS + νγ +η dt + 2 ∂S 2 ∂ν 2 ∂r2 √ √ ∂2V ∂2V ∂2V + ρ12 νγS dt + ρ13 νηS + ρ23 νγη ∂S∂ν ∂S∂r ∂ν∂r

(70)

Put-Call pariteit

We kunnen de put en call optie van eenzelfde aandeel met elkaar verbinden dankzij de Put-Call pariteit. De waarde van een call optie op expiratietijd kan beschreven worden als: C(S, T ) = max(S − K, 0)

(71)

waarbij K de strike is en T de experiratiedtijd. We zien dus, zoals verwacht, dat de call optie alleen geldt oplevert als de aandeelwaarde boven de exercise prijs ligt. Analoog vinden we voor de prijs van een put optie P (S, T ) = max(K − S, 0).

(72)

We kunnen deze formules gebruiken om put en call opties van eenzelfde aandeel met elkaar te verbinden. Stel dat we een portfolio hebben met een aandeel, een put optie en short een call optie. De call en put optie hebben dezelfde experiatijd T en strike K. We gebruiken Π voor de waarde van dit portfolio, dan geldt Π = S + P − C.

(73)

De uitbetaling van dit portfolio op expiratietijdstip is S + max(K − S, 0) − max(S − K, 0).

(74)

We onderscheiden nu 2 gevallen, S ≤ K en S ≥ K, dan vinden we:

Π(S, T ) = S + (K − S) − 0 = K als S ≤ K Π(S, T ) = S + 0 − (S − K) = K als S ≥ K We zien dus dat het portfolio altijd K waard is op tijdstip T , ongeacht het gedrag van het onderliggende aandeel. Door te verdisconteren met de rente r, vinden we de waarde van Π voor een willekeurig tijdstip t: 47

S + P − C = Ke−r(T −t) .

(75)

Relatie (75) wordt de Put-Call pariteit genoemd.

C

Upwind schema

Bij een upwindschema maken we gebruik van eenzijdige differenties om eerste orde afgeleides te benaderen, waarbij de richting van de differentie bepaald wordt door de coëfficient voor de afgeleide. Als voorbeeld behandelen we de Heston vergelijking (52). Bij centrale differentie gebruikten we voor ∂V ∂S i i Vj+1,m − Vj−1,m ∂V = ∂S 2h

met de bekende indices. Bij upwind gebruiken we een van de volgendee twee eenzijdige differenties:

∂V ∂S ∂V ∂S

= =

i i Vj+1,m − Vj,m , h i −Vi Vj,m j−1,m . h

(76) (77)

De keuze hangt af van de term voor de afgeleide. In dit geval is dit rS (merk op dat we de ∂V vergelijking bekijken in de vorm ∂V ∂t + rS ∂S + . . . = 0). Als rS > 0 maken we gebruik van (76), als rS < 0 maken we gebruik van (77). Het doel van het upwind schema is instabiliteiten uit de oplossing te halen.

D D.1 % % % %

MATLAB Appendix BlackScholesCall.m

BlackScholesCall is een functie die met behulp van de Black Scholes vergelijking de waarde van een Call optie numeriek berekent. We volgen de methode zoals gegeven in Hoofdstuk 24 in Higham, met voorwaartse differentie in de tijd.

% Input strike K, wortel van de volatiliteit sigma, rente r, expiratietijd % T, max waarde aandeel Smax, gridsize in S NS, gridsize in t Nt % Output matrix V met Black-Scholes waarde, plot van V tegen S,t function[V] = BlackScholesCall(K,sigma,r,T,Smax,NS,Nt)

48

% we hebben de substitutie tau=T-t gebruikt, waardoor we een beginwaarde % probleem hebben. % we gebruiken de relatie V(i+1)= F V(i)+p(i), met F en p als in p. 258-259 % in Higham % V is een vector van lengte NS-1 (zonder de randen S=0 en S=Smax). % We maken een matrix V aan waarin elke kolom i V(i) bevat. tau=0 geeft i=1, V is dus % een matrix van grootte NS-1 x Nt+1. % Parametersuggesties Oosterlee: sigma : 0.15 -- .3; r: .05 V = zeros(NS-1,Nt+1); % h is stapgrootte in S h = Smax/NS; % k is stapgrootte in t k = T/Nt; % Invullen beginvoorwaarde Sj=h:h:Smax-h; V(:,1) = max(Sj-K,0); % Uitrekenen F (zie p. 258 in Higham) D1 = diag([1:NS-1]); D2 = diag([1:NS-1].^2); T1 = diag(ones(NS-2,1),1) - diag(ones(NS-2,1),-1); T2 = diag(ones(NS-2,1),1) + diag(ones(NS-2,1),-1) - 2*eye(NS-1); F = (1-r*k)*eye(NS-1) + 1/2*k*sigma^2*D2*T2 + 1/2*k*r*D1*T1; % Uitrekenen P (zie p.259 in Higham) size = (NS-1,1) %P = [zeros(NS-2,1); 1/2*k*(NS-1)*(sigma^2*(NS-1)+r)*Smax]; P = [zeros(NS-2,1); 1/2*k*(NS-1)*(sigma^2*(NS-1)+r)*(Smax - K)];

for i = 1:Nt V(:,i+1) = F*V(:,i) + P; end %Toevoegen begin en eindvoorwaarde in S beginwaarde = zeros(1,Nt+1); % Eindwaarde = C(Smax,T - tau) = Smax - Ke^{-r tau) %eindwaarde = Smax*ones(1,Nt+1); %eindwaarde = Smax - K* exp(-r*(T-[0:k:T])); %eindwaarde = Smax - K* exp(-r*([0:k:T])); eindwaarde = (Smax-K)*ones(1,Nt+1); 49

V=[beginwaarde;V;eindwaarde]; %Plotten mesh([0:k:T],[0:h:Smax],V) xlabel(’T-t’),ylabel(’S’),zlabel(’Call Waarde’)

D.2 % % % %

BSExactCallPlot.m

ExactCallPlot plot de exacte Black-Sholes-waarde van een Call-optie voor verschillendeS,t. Hij gebruikt hiervoor de eerder aangemaakte functie BlackScholesExactCall en berekent deze voor een grid zoals in BlackScholesCall.

% Input strike K, wortel van de volatiliteit sigma, rente r, expiratietijd % T, eindwaarde aandeel Smax, gridsize in S NS, gridsize in t Nt % Output matrix V met Black-Scholes waarde, plot van V tegen S,t function[V] = BSExactCallPlot(K,sigma,r,T,Smax,NS,Nt) V = zeros(NS+1,Nt+1); for i=1:NS+1 for j=1:Nt+1 V(i,j) = BlackScholesExactCall((i-1)/NS*Smax,K,r,sigma, (j-1)/Nt * T); end end mesh([0:T/Nt:T],[0:Smax/NS:Smax],V) xlabel(’T-t’),ylabel(’S’),zlabel(’Call Waarde’)

D.3

BlackScholesExactCall.m

%BlackScholesExactCall rekent de exacte waarde uit van een Call optie. De %code komt rechtstreeks uit het boek van Higham (CH.08 p. 84) %Input: asset price S, strike K, rente r, volatiliteit sigma, tijd tot %expiratiedatum t %Output: Call waarde C function[C] = BlackScholesExactCall(S,K,r,sigma,tau) if tau>0 d1 = (log(S/K) + (r+.5*sigma^2)*(tau))/(sigma*sqrt(tau)); d2 = d1 - sigma*sqrt(tau); N1 = 0.5*(1+erf(d1/sqrt(2))); N2 = 0.5*(1+erf(d2/sqrt(2))); C = S*N1-K*exp(-r*(tau))*N2; 50

else C=max(S-K,0); end

D.4 % % % %

TestBlackScholes.m

TestBlackScholes plot de exacte Black-Sholes-waarde van een Call-optie voor verschillende S,t. Hij gebruikt hiervoor de eerder aangemaakte functie BlackScholesExactCall en berekent deze voor een grid zoals in BlackScholesCall.

% Input strike K, wortel van de volatiliteit sigma, rente r, expiratietijd % T, eindwaarde aandeel Smax, gridsize in S NS, gridsize in t Nt % Output matrix V met Black-Scholes waarde, plot van V tegen S,t function[Vnum,Vexact,VerschilMatrix,MaxVerschil,GemVerschil,stabi,tijd] = TestBlackScholes(NS,Nt) tic; close all % Parameters K = 100; sigma = .2; r = 0.05; T = 1; Smax = 200; % Numerieke oplossing figure(’name’,’Numerieke waarde’) Vnum = BlackScholesCall(K,sigma,r,T,Smax,NS,Nt); % Exacte oplossing Vexact = zeros(NS+1,Nt+1); for i=1:NS+1 for j=1:Nt+1 Vexact(i,j) = BlackScholesExactCall((i-1)/NS*Smax,K,r,sigma, (j-1)/Nt * T); end end % Plotten exacte oplossing figure(’name’,’Exacte oplossing’) mesh(0:T/Nt:T,0:Smax/NS:Smax,Vexact) xlabel(’T-t’),ylabel(’S’),zlabel(’Call Waarde’) % Verschil tussen numerieke en exacte oplossing VerschilMatrix = Vnum - Vexact; 51

VerschilMatrixAbs = abs(VerschilMatrix); %Absolute waardes % Plotten VerschilMatrixAbs figure(’name’,’Absolute fout numerieke schema’) mesh(0:T/Nt:T,0:Smax/NS:Smax,VerschilMatrixAbs) xlabel(’T-t’),ylabel(’S’),zlabel(’Fout’) % Plotten VerschilMatrix figure(’name’,’Fout numerieke schema’) mesh(0:T/Nt:T,0:Smax/NS:Smax,VerschilMatrix) xlabel(’T-t’),ylabel(’S’),zlabel(’Fout’) % Maximale verschil MaxVerschil = max(max(VerschilMatrixAbs)); % Gemiddeld absoluut verschil GemVerschil = sum(sum(VerschilMatrixAbs))/(size(VerschilMatrixAbs,1)*size(VerschilMatrixAbs,2)); % Plotten fout op t=0 figure(’name’,’Fout op t=0’) plot(0:Smax/NS:Smax,VerschilMatrixAbs(:,size(VerschilMatrixAbs,2))) xlabel(’S’),ylabel(’Fout’) stabi = (T/Nt)/(Smax/NS); tijd = toc; end

D.5

HestonCall.m

% HestonCall past een numeriek schema toe op het Heston model met voorwaartse % differenties in de tijd en centrale differenties in v en S. De code wordt % beschreven in het BSc project van Pim Stuurman % % % % % %

Output: Vmatrix: matrix die de optiewaarde op tau=0 (t=T) bevat. De matrix heeft grootte (NS+1)x(Nv+1); looptijd: runtime in seconden; testwaarde: optiewaarde op tau=T (t=0), v=0.04 en S=100; 3d plot van optiewaarde tegen v en S op eindtijd tau=T (t=0).

% Nummering (v=m,S=j) -> m+1 + j(Nv+1), waarbij Nv het grootste element van % v is. function[Vmatrix,looptijd,stabi] = 52

HestonCall(Smax,vmax,NS,Nt,Nv,testcase,randvoorwaarde) tic; % Gegevens voor alle testcases K = 100; r = 0; rho = -0.9; vstreep = 0.04; % Gegevens voor specifieke testcases if testcase == 1 gamma = 0.5; kappa = 5; T=1; elseif testcase == 2 gamma = 0.5; kappa = 0.5; T=1; elseif testcase == 3 gamma = 1; kappa = 0.5; T=10; end % V is een vector van grootte (Nv+1)x(NS+1). V wordt in de uiteindelijke % for-loop recursief berekend. Aan het end % Aanmaken V. V= zeros((Nv+1)*(NS+1),1); % l h k

Stapgroottes in v,S en t respectievelijk. = vmax/Nv; = Smax/NS; = T/Nt;

% % % %

Aanmaken D_i (de randen van het grid worden later uit de matrices gehaald) We maken eerst een matrix aan met de goede elementen, die ordenen we verder volgens horizontale nummering in de v-richting.

% D1 = mj^2 elementen1 = [0:Nv]’*[0:NS].^2; D1 = diag(elementen1(:)); % D2 = mj elementen2 = [0:Nv]’*[0:NS]; D2 = diag(elementen2(:)); % D3 = m elementen3=[0:Nv]’*ones(NS+1,1)’; 53

D3 = diag(elementen3(:)); % D4 = j elementen4=[0:NS]’*ones(Nv+1,1)’; D4 = diag(elementen4(:)); % D5 = m elementen5=[0:Nv]’*ones(NS+1,1)’; D5 = diag(elementen5(:)); % Aanmaken T_i (de randen van het grid worden later uit de matrices % gehaald) T1= diag(ones((Nv+1)*NS,1),(Nv+1)) -2*eye((Nv+1)*(NS+1)) + diag(ones((Nv+1)*NS,1),-(Nv+1)); T2= diag(ones((Nv+1)*NS-1,1),(Nv+2)) - diag(ones((Nv+1)*NS+1,1),-Nv) diag(ones((Nv+1)*NS+1,1),Nv) + diag(ones((Nv+1)*NS-1,1),-(Nv+2)); T3= diag(ones((Nv+1)*(NS+1)-1,1),1) - 2*eye((Nv+1)*(NS+1)) + diag(ones((Nv+1)*(NS+1)-1,1),-1); T4= diag(ones((Nv+1)*NS,1),(Nv+1)) - diag(ones((Nv+1)*NS,1),-(Nv+1)); T5= diag(ones((Nv+1)*(NS+1)-1,1),1) - diag(ones((Nv+1)*(NS+1)-1,1),-1); T6= T5; % F opstellen zonder randvoorwaarden F = 1/2*k*l*D1*T1 + 1/4*rho*gamma*k*D2*T2 + 1/(2*l)*gamma^2*k*D3*T3 + 1/2*r*k*D4*T4 - (r*k-1)*eye((Nv+1)*(NS+1)) - 1/2*k*kappa*D5*T5 + 1/(2*l)*k*vstreep*kappa*T6; % Vervangen elementen op rand v=vmax. Dit zijn de laatste Nv+1 rijen in de % matrix. Omdat het een Dirichlet voorwaarde is, kunnen we de elementen % uit D_i en T_i overal vervangen door de identiteit. vervanging_vmax = eye((Nv+1)*(NS+1)); F(NS*(Nv+1)+1:(NS+1)*(Nv+1),:) = vervanging_vmax(NS*(Nv+1)+1:(NS+1)*(Nv+1),:); % Vervangen elementen op rand S=0 vervanging_S0 = eye((Nv+1)*(NS+1)); F(1:Nv+1:(Nv+1)*(NS+1),:) = vervanging_S0([1:Nv+1:(Nv+1)*(NS+1)],:); % Vervangen elementen op rand S=Smax vervanging_Smax = eye((Nv+1)*(NS+1)); F(Nv+1:Nv+1:(Nv+1)*(NS+1),:) = vervanging_Smax([Nv+1:Nv+1:(Nv+1)*(NS+1)],:); % % % % % % %

Vervangen elementen op rand v=0. Verschillende keuzes voor verschillende randvoorwaardes. Overzicht: -1 = heston -2 = alternatief 1 -3 = alternatief 2 enz. Voor beschrijving randvoorwaarden zie verslag.

% Om eerste regels van eye((Nv+1)*(NS+1)) te kunnen kiezen, maken we deze eerst aan 54

I=eye((Nv+1)*(NS+1)); if randvoorwaarde == -1 vervanging_v0 = 1/2*k*r*D4(1:Nv+1,:)*T4 + 1/(2*l)*k*vstreep*kappa*T5(1:Nv+1,:)+(1-r*k)*I(1:Nv+1,:); F(1:Nv+1,:) = vervanging_v0(1:Nv+1,:); else disp(’FOUT: VOER GELDIGE RANDVOORWAARDE KEUZE IN!’); end % Invullen beginwaarde tau=0. Sj=0:h:Smax; beginvoorwaarde = max(Sj-K,0); for i=0:NS % V([i*(Nv+1)+1:(i+1)*(Nv+1)],1) = beginvoorwaarde(i+1); V(i*(Nv+1)+1:(i+1)*(Nv+1)) = beginvoorwaarde(i+1); end % Uitvoeren schema V= F^(Nt)*V; Vmatrix = vec2mat(V,Nv+1); %Plotten mesh(0:l:vmax,0:h:Smax,Vmatrix) xlabel(’v’),ylabel(’S’),zlabel(’Call Waarde’) % Stabiliteitsindicator stabi=k/(l^2+h^2); looptijd=toc; end % % % % % % % % % % % % % % % % %

% % % % % % % % % % % % % % % % %

% De test wordt geevalueerd op het punt v=0.04, S=100, T=1. We zoeken die % waarde op. m en j zijn de indices in v en S, respectievelijk. m=(0.04*Nv)/vmax + 1; j=(100*NS)/Smax + 1; % Als m of j geen geheel getal is, moeten we afronden. Dan komt er een % extra fout tov de exacte oplossing, omdat we niet meer in het goede punt % zoeken. Er wordt hier melding van gegeven.

if rem(m,1) ~= 0 m=round((0.04*Nv)/vmax + 1); disp(’Er moet afgerond worden in volatiliteitsindex, we introduceren extra fout!’) elseif rem(j,1) ~=0 j=round((100*NS)/Smax + 1); disp(’Er moet afgerond worden in Stockindex, we introduceren extra fout!’); end 55

% % % %

% % %In Matlab gaat de indicatie andersom dan in verslag. Daarom staat er % %Vmatrix(j,m) ipv (m,j). % testwaarde = Vmatrix(j,m);

D.6 % % % %

HestonCallIte.m

HestonCallIte past het numerieke schema toe om de waarde van een optie in het Heston model numeriek te berekenen. Het is een iteratieve versie van de file HestonCall.m. Verdere toelichting is te vinden in die file en in het verslag.

function[Vmatrix,looptijd,stabi] = HestonCallIte(Smax,vmax,NS,Nt,Nv,testcase,randvoorwaarde) tic; % Gegevens voor alle testcases K = 100; r = 0; rho = -0.9; vstreep = 0.04; % Gegevens voor specifieke testcases if testcase == 1 gamma = 0.5; kappa = 5; T=1; elseif testcase == 2 gamma = 0.5; kappa = 0.5; T=1; elseif testcase == 3 gamma = 1; kappa = 0.5; T=10; end % V is een vector van grootte (Nv+1)x(NS+1). V wordt in de uiteindelijke % for-loop recursief berekend. % Aanmaken V. V= zeros((Nv+1)*(NS+1),1); % Stapgroottes in v,S en t respectievelijk. l = vmax/Nv; h = Smax/NS; 56

k = T/Nt; % Invullen beginwaarde tau=0. Sj=0:h:Smax; % % % % % %

%Smoothing rond het punt S=K Smoothingindex = K/h+1; %locatie van het knikpunt (S(index) is laatste 0) SmoothN = 2; %aantal punten om knik gladder te maken for i=0:SmoothN-1 Sj(Smoothingindex-i) = (SmoothN-i/(SmoothN+1))*Sj(Smoothingindex+1); end

if CallPut ==1 %Call beginvoorwaarde = max(Sj-K,0); else %Put beginvoorwaarde = max(K-Sj,0); end for i=0:NS V(i*(Nv+1)+1:(i+1)*(Nv+1)) = beginvoorwaarde(i+1); end % Invullen randvoorwaardes (v=0 is verwerkt in for-loop) % Rand v=vmax %for j=0:NS % n=Nv+1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). % V(n)= max(Sj(j+1)-K,0); %V(n)= Sj(j+1); %end % Rand S=0 if CallPut ==1 % Call for m=0:Nv n=m+1; %indicator in horizontale nummering. V(n)= 0; end end %Put wordt behandeld in de for-loop % Rand S=Smax for m=0:Nv n=m+1+NS*(Nv+1); %indicator in horizontale nummering. %V(n)=V(n-(Nv+1)) + h; %dV/dS=1 if CallPut == 1 % Call V(n)= Smax-K; %V(n)=Smax; else %Put V(n)=0; end end %For-loop voor eindige differenties 57

for i = 1:Nt % Rand v=0 % keuzes worden beschreven in verslag % randvoorwaarde ==-1 is de max(S-K,0), deze is standaard al verwerkt. if randvoorwaarde==-2 for j=1:NS-1 n=1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). V(n)= (k*r*j)/2*(V(n+Nv+1)-V(n-(Nv+1))) (k*kappa*vstreep)/(2*l)*(V(n+1)-V(n-1))+(1-r*k)*V(n); end elseif randvoorwaarde==-3 %Voorwaartse differenties for j=1:NS-1 n=1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). V(n)= (k*r*j)*(V(n+Nv+1)-V(n)) (k*kappa*vstreep)/l*(V(n+1)-V(n))+(1-r*k)*V(n)+ (gamma^2*m)/2*(V(n+1) - V(n)); end elseif randvoorwaarde==-4 for j=1:NS-1 n=1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). V(n)= (k*r*j)/2*(V(n+Nv+1)-V(n-(Nv+1))) (k*kappa*vstreep)/(2*l)*(V(n+1)-V(n-1))+(1-r*k)*V(n)+ (gamma^2*m)/(2*l)*(V(n+1) - 2*V(n)+V(n-1))+ (m*l*j^2)/2*(V(n+(Nv+1))-2*V(n)+V(n-(Nv+1))); end elseif randvoorwaarde==-5 for j=1:NS-1 n=1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). % randvoorwaarde uit boek %Upwind if kappa*vstreep>0 V(n)=k*r*j*(V(n+Nv+1)-V(n))+ k*kappa*vstreep*(V(n+1)-V(n))/l +(1-r*k)*V(n) else V(n)=k*r*j*(V(n+Nv+1)-V(n))+ k*kappa*vstreep*(V(n)-V(n-1))/l +(1-r*k)*V(n) end end end % Rand S=0 voor Put if CallPut == 2 %Put for m=0:Nv n=m+1; %indicator in horizontale nummering. V(n)= K*exp(-r*k*i); end end %interne elementen for m = 1:Nv-1 for j=1:NS-1 58

n=m+1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). V(n) = (k*m*l*j^2)/2 * (V(n+Nv+1) - 2*V(n) + V(n-(Nv+1))) + (rho*gamma*k*m*j)/4*(V(n+Nv+2)-V(n-(Nv+1)+1)-V(n+(Nv+1)-1)+V(n-Nv-2)) + (gamma^2*k*m)/(2*l)*(V(n+1)-2*V(n)+V(n-1))+ (r*k*j)/2*(V(n+Nv+1)-V(n-(Nv+1)))- (r*k-1)*V(n) (k*m*kappa)/2*(V(n+1)-V(n-1))+ (k*vstreep*kappa)/(2*l)*(V(n+1)-V(n-1)); % Upwind in dV\dv: % if ((vstreep*kappa)-(l*m*kappa))<0 % V(n) = (k*m*l*j^2)/2 * (V(n+Nv+1) - 2*V(n) + V(n-(Nv+1))) + % (rho*gamma*k*m*j)/4*(V(n+Nv+2)-V(n-(Nv+1)+1)-V(n+(Nv+1)-1)+V(n-Nv-2)) + % (gamma^2*k*m)/(2*l)*(V(n+1)-2*V(n)+V(n-1))+ (r*k*j)/2*(V(n+Nv+1)-V(n-(Nv+1)))% (r*k-1)*V(n) - (k*m*kappa)*(V(n+1)-V(n))+ (k*vstreep*kappa)/(l)*(V(n+1)-V(n)); % else % V(n) = (k*m*l*j^2)/2 * (V(n+Nv+1) - 2*V(n) + V(n-(Nv+1))) + % (rho*gamma*k*m*j)/4*(V(n+Nv+2)-V(n-(Nv+1)+1)-V(n+(Nv+1)-1)+V(n-Nv-2)) + % (gamma^2*k*m)/(2*l)*(V(n+1)-2*V(n)+V(n-1))+ % (r*k*j)/2*(V(n+Nv+1)-V(n-(Nv+1)))- (r*k-1)*V(n) - (k*m*kappa)*(V(n)-V(n-1))+ % (k*vstreep*kappa)/(l)*(V(n)-V(n-1)); % end end end % rand v=vmax (dV/dv =0) met voorwaartse differentie. for j=0:NS n=Nv+1+j*(Nv+1); %indicator voor horizontale nummering: V(v=m,S=j)=V(n). V(n)= V(n-1); end

% % % %

% rand S=smax (dV/ds=0) met voorwaartse differentie. for m=0:Nv n=m+1+NS*(Nv+1); %indicator in horizontale nummering. V(n)= V(n-(Nv+1)); end

end Vmatrix = vec2mat(V,Nv+1); %Plotten mesh(0:l:vmax,0:h:Smax,Vmatrix) xlabel(’v’),ylabel(’S’),zlabel(’Call Waarde’) % Stabiliteitsindicator stabi=k/(l^2+h^2); looptijd=toc; end

59

D.7

% % % %

HestonExactCall.m

HestonExactCall bepaalt de exacte oplossing van het Heston model. De code komt uit ’Finite Difference Schemes for Heston Model’ van S. Lin. De parameters zijn niet exact gelijk aan die in onze notatie, ze kunnen uiteraard in elkaar omgeschreven worden.

function R3 = HestonExactCall(S,v,testcase,trunc) % Gegevens voor alle testcases K = 100; r = 0; rh = -0.9; et = 0.04; ld=0; %ld is lambda, deze ontstaat door de scheiding van variabelen. % Gelijk aan theta in mijn BSc verslag. greek = 1; %Met deze code kunnen ook andere greeks gekozen worden, % maar we zijn alleen geinteresseerd in de prijs. % Gegevens voor specifieke testcases if testcase == 1 sm = 0.5; kp = 5; t=1; elseif testcase == 2 sm = 0.5; kp = 0.5; t=1; elseif testcase == 3 sm = 1; kp = 0.5; t=10; end % Oorspronkelijke code: if nargin < 11, trunc = 100; greek=1; end I = sqrt(-1); x = log(S)+r*t; switch greek case 1 % price inte = @(w) (exp(-I*w*x).*(K.^(1+I*w)).*Ffun(t,w,v,kp,ld,et,sm,rh)./(I*w-w.^2)); case 2 % delta inte = @(w) (-I*w./S.*exp(-I*w*x).*(K.^(1+I*w)).*Ffun(t,w,v,kp,ld,et,sm,rh)./(I*w-w.^2));

60

case 3 % gama inte = @(w) (-(w.^2)./(S^2).*exp(-I*w*x).*(K.^(1+I*w)).*Ffun(t,w,v,kp,ld,et,sm,rh)./(I*w-w.^2)); case 4 % vega inte = @(w) (exp(-I*w*x).*(K.^(1+I*w)).*Ffun(t,w,v,kp,ld,et,sm,rh,1)./(I*w-w.^2)); end R3t = exp(-r*t)*quadgk(inte,-trunc+2i,trunc+2i)/(2*pi); R3 = real(R3t); end function R1 = Ffun(t,w,v,kp,ld,et,sm,rh,Indi) if nargin<9 Indi = 0; end I d g D C

= = = = =

sqrt(-1); sqrt((rh*sm*I*w+kp+ld).^2+(w.^2-I*w).*(sm^2)); (kp+ld+rh*sm*I*w+d)./(kp+ld+rh*sm*I*w-d); ((kp+ld+rh*sm*I*w+d)/(sm^2)).*(1-exp(d*t))./(1-g.*exp(d*t)); kp*et*((kp+ld+rh*sm*I*w+d)*t-2*log((1-g.*exp(d*t))./(1-g)))/(sm^2);

if Indi ==0 R1 = exp(C+D.*v); else R1 =D.*exp(C+D.*v); end end

D.8

TestHeston.m

% TestHeston vergelijkt de numerieke resultaten van HestonCall.m met de % exacte oplossing uit HestonExactCall.m % % % %

Ideeen: plots van beide gevallen, verschilplot, maximale verschil (en plaats van verschil), s nachts laten runnen over verschillende verhoudingen van NS,Nv en Nt bekijk verschillende testcases.

function [Vnum,Vexact,VerschilMatrix,MaxVerschil,GemVerschil,stabi,tijd] = TestHeston(Smax,vmax,NS,Nt,Nv,testcase,randvoorwaarde) tic; close all 61

% Aanroepen HestonCall. Hierbij wordt een vector gegenereerd met de waarden % op t=0 (tau=T) figure(’name’,’Numerieke waarde’) [Vnum,~,stabi] = HestonCallIte(Smax,vmax,NS,Nt,Nv,testcase,randvoorwaarde); % Aanmaken van een matrix met de exacte oplossingen. De vector moet % dezelfde punten analyseren als Vmatrix. % Aanmaken lege matrix Vexact = zeros((NS+1),(Nv+1)); % l h v S

Aanmaken vectoren S en v waarin de te berekenen punten zitten = vmax/Nv; = Smax/NS; = 0:l:vmax; = 0:h:Smax;

% Invullen exacte oplossing: Vexact(1,:)=0; %Voor S=0 vullen we handmatig in, daarvoor is geen exacte oplossing for i=1:Nv+1 for j=2:NS+1 Vexact(j,i) = HestonExactCall(S(j),v(i),testcase,100); end end % Plotten Vexact figure(’name’,’Exacte Waarde’) mesh(0:l:vmax,0:h:Smax,Vexact) xlabel(’v’),ylabel(’S’),zlabel(’Call Waarde’) % Verschil tussen numerieke en exacte oplossing VerschilMatrix = Vnum - Vexact; VerschilMatrixAbs = abs(VerschilMatrix); %Absolute waardes % Plotten VerschilMatrixAbs figure(’name’,’Absolute waarde fout numerieke schema’) mesh(0:l:vmax,0:h:Smax,VerschilMatrixAbs) xlabel(’v’),ylabel(’S’),zlabel(’Fout’) % Plotten VerschilMatrixAbs figure(’name’,’Fout numerieke schema’) mesh(0:l:vmax,0:h:Smax,VerschilMatrix) xlabel(’v’),ylabel(’S’),zlabel(’Fout’) % Maximale verschil MaxVerschil = max(max(VerschilMatrixAbs)); 62

% Gemiddeld verschil GemVerschil = sum(sum(VerschilMatrixAbs))/(size(VerschilMatrixAbs,1)*size(VerschilMatrixAbs,2)); tijd=toc; end

63

Technische Universiteit Delft Faculteit Elektrotechniek, Wiskunde en Informatica Delft Institute of Applied Mathematics

Recommend Documents