Equacions en derivades parcials a les finances

Abstract Options trading is a fundamental part of finance. The Black-Scholes partial differential equation models the price of these contracts. For many options, there are well-known solutions to the equation but numerical problems arise with American style options and other type of options due to the lack of analytical solutions. This project aims to provide a mathematical model to solve numerically this type of pricing problems. Keywords finance, options, Black-Scholes, numerical PDE, numerical methods 1 ´ Index 1 Introducci´o 4 1.1 Opcions ............................................ 4 1.1.1 Altres tipus d’opcions i opcions ex`otiques ...................... 5 1.2 Tipus d’inter`es lliure de risc i arbitratge ........................... 7 2 L’equaci´o de Black-Scholes 9 2.1 Formulaci´o ........................................... 9 2.2 Opcions americanes ...................................... 12 2.3 Opcions asi`atiques ....................................... 14 3 M`etodes num`erics per l’equaci´o de Black-Scholes 16 3.1 M`etodes de difer`encies finites ................................. 16 3.1.1 Aproximaci´o de les derivades de V(s,t)....................... 16 3.2 Opcions europees ....................................... 17 3.2.1 M`etode Expl´ıcit .................................... 17 3.2.2 M`etode Impl´ıcit .................................... 19 3.2.3 M`etode Crank-Nicolson ................................ 21 3.2.4 Resoluci´o de sistemes lineals ............................. 22 3.3 M`etodes de Montecarlo .................................... 23 3.3.1 T`ecnica reducci´o vari`ancia .............................. 24 3.3.2 Exemple: Opcions europees .............................. 25 4 Calcul d’opcions americanes i asi`atiques 26 4.1 Opcions americanes ...................................... 26 4.2 Opcions asi`atiques ....................................... 27 4.3 Comparaci´o dels preus entre diferents tipus d’opcions i efecte dels diferents par`ametres en el preu ............................................. 28 5 Bibliografia 30 A Codis M`etode Diferencies Finites 31 A.1 M`etode Expl´ıcit ........................................ 31 A.2 M`etode Impl´ıcit ........................................ 32 A.3 M`etode Crank-Nicolson .................................... 33 A.4 M`etode opcions americanes .................................. 35 A.5 M`etode de sobre-relaxaci´o successiva per a sistemes lineals ................. 35 2 B Codi M`etode de Montecarlo 36 B.1 M`etode opcions europees ................................... 36 B.2 M`etode opcions asi`atiques .................................. 36 3 Edps a les finances 1. Introducci´o L’objectiu d’aquest treball de fi de grau ´es resoldre per diferents m`etodes num`erics l’equaci´o de BlackScholes que regula el preu de les opcions als mercats financers. Primerament introduirem els elements del m´on financer que intervenen en l’equaci´o de Black-Scholes, que seguidament derivarem usant regles pr`opies del m´on de les finances (per exemple el concepte d’arbitratge). Despr`es resoldrem l’equaci´o amb diversos m`etodes num`erics per a poder finalment fer una comparaci´o dels preus de diversos tipus d’opcions. En aquest cap´ıtol introdu¨ım el concepte d’opcions europees, americanes i asi`atiques que m´es tard seran l’objecte d’estudi i en modelarem el seu preu a trav´es de resoldre num`ericament l’equaci´o de Black-Scholes. Tamb´e introduirem conceptes b`asics per al seu estudi com poden ser el tipus d’inter`es lliure de risc o l’arbitratge. 1.1 Opcions Una opci´o ´es un producte financer derivat, ´es a dir, que dep`en d’un altre producte financer subjacent, per exemple les accions d’una empresa o el preu d’una mat`eria prima. Definici´o 1.1. Una opci´o ´es un contracte entre dues parts en qu`e el comprador (posici´o llarga) adquireix al venedor (posici´o curta) el dret de comprar (o vendre) el producte subjacent a un cert preu fixat Ken una data fixada T.Ks’anomena el preu d’exercici (strike price) i Ts’anomena el temps de venciment (expiration date). Observaci´o 1.2.El comprador de l’opci´o t´e el dret, per`o no l’obligaci´o de comprar (o vendre) el subjacent. Observaci´o 1.3.Si el comprador adquireix el dret a comprar el subjacent estarem parlant d’una opci´o call, en cas de venda parlarem d’una opci´o put. Distingim ara entre opcions europees i americanes. Definici´o 1.4. Una opci´o europea ´es una opci´o que nom´es pot ser executada a data de venciment. Exemple 1.5. Un exemple d’un call europeu ´es un contracte entre AiBon Aadquireix a Bel dret a comprar 100 accions d’una companyia (subjacent) a un preu de K= 2 per acci´o (strike price) d’aqu´ı un any (expiration date), per tant T= 1 any. Passat un any At´e l’opci´o de comprar el paquet d’accions al preu acordat, per`o nom´es ho far`a si les accions valen en el mercat m´es del preu acordat, de no ser aix´ı les podria comprar a mercat a menor preu i el contracte no val res. Definici´o 1.6. Una opci´o americana ´es una opci´o que pot ser executada en qualsevol moment des que es firma el contracte fins a temps de venciment. Exemple 1.7. Si ens posem en el mateix cas que en l’exemple anterior la difer`encia amb una opci´o europea est`a en el fet que Apot executar la opci´o (comprar el paquet d’accions al preu acordat) a temps t<T, si veu que en cert moment abans de la data de venciment obt´e benefici. Molts cops es pensa que les opcions nom´es serveixen com a eina especuladora. Les opcions tamb´e es poden usar per a cobrir el risc de certes operacions (aix´ı com ho faria una asseguran¸ca). 4 Exemple 1.8. Suposem que l’empresa A ven blat (subjacent) al mercat. Si el mercat actual est`a a 1/kg i l’empresa A espera vendre el blat d’aqu´ı un any pot adquirir un put europeu a K= 0.95/kg. A temps de venciment, si el mercat est`a per sobre del preu d’exercici A no executar`a el put i vendr`a el blat al mercat, nom´es l’executar`a si el preu del blat esta per sota del strike price. D’aquesta manera A es protegeix d’una eventual baixada del preu del blat. Tractem ara la liquidaci´o d’una opci´o (anomenat payoff ). Sigui STel preu del subjacent a temps de venciment i Kel seu preu d’exercici. Si s’executa una opci´o de compra s’obt´e ST−K, de la mateixa manera, en executar una opci´o de venda s’obt´e K−ST. Clarament les opcions nom´es s’executen si la seva liquidaci´o es positiva. D’aquesta manera, i sent sel preu del subjacent, denotem per φ(s) la liquidaci´o d’una opci´o: φ(s) =    max(ST−K, 0) (compra) max(K−ST, 0) (venda) (a) Liquidaci´o opci´o de compra (b) Liquidaci´o opci´o de venda Figura 1: Liquidacions per a opcions amb K= 100 1.1.1 Altres tipus d’opcions i opcions ex`otiques A part de les opcions que hem explicat es poden crear altres tipus d’opcions m´es complexes que les explicades anteriorment. Per exemple, suposem que comprem una call a preu d’exercici 100 i venem una call a preu d’exercici 120. Aleshores, acabem de crear una cartera d’opcions anomenat bull spread. En la seg¨uent figura mostrem el payoff de la nostra cartera. La f´ormula general del payoff d’un bull spread amb una call comprada de preu d’exercici K1i una venuda de preu d’exercici K2(on K2>K1) ´es φ(s) = 1 K2−K1 (max(s−K1, 0) −max(s−K2, 0)) De la mateixa manera que hem constru¨ıt aquesta opci´o composta es pot arribar a construir una cartera d’opcions amb la funci´o de payoff desitjada. ´ Es un proc´es que no explicarem de manera detallada per`o que involucra una cartera d’opcions amb una idea similar a l’emprada per a crear les anomenades opcions de barrera. De totes formes modificar la funci´o de liquidaci´o no ´es l’´unica manera de construir opcions m´es complexes. Una de les m´es conegudes s´on les opcions de barrera que consisteixen en opcions de compra i/o venda 5 Edps a les finances Figura 2: Liquidaci´o bull spread est`andards amb l’afegit que si el preu del subjacent assoleix un cert preu (preu de barrera) l’opci´o perd tot el valor. D’aquesta manera s’introdueix un nou concepte, el de la depend`encia de cam´ı, ja que es dona el cas que dos camins (entesos com l’evoluci´o del valor del subjacent respecte el temps) que a temps d’exercici tinguin el mateix valor poden no tenir mateixa liquidaci´o. Com podem veure en la figura (13) a temps de venciment el preu del subjacent est`a per sobre el preu Figura 3: Evoluci´o del preu d’un subjacent en una call de strike K = 120 amb barrera en S= 90. d’exercici. De totes maneres l’opci´o no val res ja que abans el preu ha passat per sota el preu de barrera i per tant l’opci´o perd tot el seu valor en aquell moment. Un altre tipus d’opcions amb depend`encia de cam´ı s´on les opcions asi`atiques. En una opci´o asi`atica la liquidaci´o dep`en de la mitjana (aritm`etica o geom`etrica) del preu del subjacent des de l’inici del contracte fins a la data de venciment en comptes de dependre del preu a venciment. Aquesta mitjana pot ser calculada de forma cont´ınua o discreta. A continuaci´o enumerem les liquidacions de les diferents opcions asi`atiques per opcions de compra (per les opcions de venda s´on molt similars) amb temps de venciment T i on el preu del subjacent durant el contracte ´es St: 1. Mitjana aritm`etica cont´ınua: 1 TRT 0Stdt −K 6 Mentre el valor de l’opci´o sigui superior al del payoff el preu de l’opci´o satisfar`a l’equaci´o de Black-Scholes, tenint en compte que el valor de l’opci´o no dep`en del temps obtenim: 1 2σ2s2∂2V ∂2s2+rs ∂V ∂s−rV = 0 si resolem aquesta equaci´o diferencial de segon ordre la seva soluci´o general ´es: V(s) = As +Bs−2r/σ2(17) Clarament A= 0 ja que si fem tendir el preu del subjacent a infinit el valor de l’opci´o (al ser una put) ha de tendir a 0. Tenim doncs V(s) = Bs−2r/σ2i suposem que existeix un valor ˆstal que exercim l’opci´o (´es a dir, tant bon punt el subjacent arribi a valdre ˆsexercirem el contracte). Tenim doncs a s= ˆsel seg¨uent: V(ˆs) = Bˆs−2r/σ2=K−ˆs De la segona igualtat podem a¨ıllar Bi si ho subtituim a (17) obtenim V(s) = (K−ˆs)s ˆs−2r/σ2 (18) Sobre la tria de ˆs, hem de triar ˆsde forma que maximitzi el valor del contracte en qualsevol instant de temps abans del temps d’exercici. Aix`o es fa ja que qui t´e l’opci´o de venda en vol maximitzar el seu valor. Derivant (18) respecte ˆsi igualant a 0 s’obt´e ˆs=K 1 + σ2/2r de forma que ja tenim V(s) completament caracteritzada. A partir d’aix`o es pot veure que la derivada del valor del contracte i del payoff es la mateixa a s= ˆs. Aix`o ens fa veure la necessitat que ∆ = ∂V ∂s sigui cont´ınua. Aquesta propietat no nom´es s’aplica en el cas de la put perp`etua sin´o en totes les opcions americanes. Finalment, les condicions (13),(14),(15)i(16) creen el problema a resoldre per a obtenir el valor de les opcions americanes. Aquest problema entra dintre de la categoria dels anomenats problemes de frontera lliure. Expliquem ara perqu`e diem que es tracta d’un problema de frontera lliure. A cada temps texisteix un valor de sque fa de frontera entre la regi´o on ´es `optim exercir el contracte i en la que no. Aquesta frontera la podem notar per Sf(t). Si Sf(t)<K(cas d’una put) aleshores el pendent del payoff al punt Sf(t) ´es −1 i un argument d’arbitratge mostra que ∂V ∂sha de ser −1 tamb´e (continu¨ıtat de la derivada). Suposem que ∂V ∂s<−1, cosa que es correspon al cas (b) de la figura (6). Aleshores V(s,t) passa per sota del payoff i aix`o contradiu (14). En el cas ∂V ∂s>−1, cas (a) de la figura (6) existeix un argument a partir de l’estrat`egia de qui t´e l’opci´o per a maximitzar-ne el seu valor que indica que no est`a ben valorada. Sf(t) ´es de fet el valor que simult`aniament maximitza el guany de la persona que t´e el contracte i n’evita la possibilitat d’arbitratge. Diem que es tracta d’un problema de frontera lliure ja que la frontera Sf(t) evoluciona amb el temps i ´es un valor desconegut a priori. En la seg¨uent llista resumim totes les condicions a resoldre per a posar preu a les opcions americanes. 13 Edps a les finances Figura 6: En aquesta figura podem observar el que passa quan (a) s’exercita l’opci´o massa tard o (b) massa d’hora. •∂V ∂t+1 2σ2s2∂2V ∂2s2+rs ∂V ∂s−rV ≤0 •∆ = ∂V ∂scont´ınua. •V(s,t)≥φ(s,t) •V(s,T) = φ(s,T) com a condici´o inicial. 2.3 Opcions asi`atiques En aquesta secci´o tractem com modelar les opcions asi`atiques usant Black-Scholes. Recordem que les opcions asi`atiques empren una funci´o de liquidaci´o que no nom´es dep`en del preu del subjacent en el moment del exercici sin´o del cam´ı que pren el preu del subjacent durant la durada de tot el contracte; sent normalment la mitjana aritm`etica (cont´ınua o discreta). En aquesta secci´o tractarem sobre una opci´o amb mitjana aritm`etica cont´ınua. Diem doncs que la liquidaci´o dep`en d’una integral d’una funci´o sobre el preu del subjacent (St) des de t= 0 fins a data de venciment t=T. En el nostre cas I(T) = ZT 0 f(s,τ)dτ=ZT 0 Sτdτ I per tant que el payoff de l’opci´o ´es una funci´o φ(s,I) a temps t=T. En el nostre cas φ(s,I) = 1 TI(T) Abans d’arribar a temps de venciment tenim informaci´o parcial sobre el preu del subjacent a temps de venciment (com m´es alt estigui Stper tabans del venciment podem esperar que STsigui m´es alt). De la mateixa manera a temps ttamb´e tenim informaci´o parcial sobre el valor de Ia temps de venciment amb 14 el valor de la integral fins a temps t: I(t) = Zt 0 f(s,τ)dτ=Zt 0 Sτdτ D’aquesta manera ens podem imaginar que el valor de l’opci´o asi`atica a temps tno nom´es dep`en de Si t, sin´o tamb´e de I.Iser`a una nova variable, anomenada variable d’estat. De la mateixa manera que el preu del subjacent segueix una equaci´o diferencia estoc`astica (com en (2)), Itamb´e ho fa. Emprant el lema de Itˆo tenim que dI =f(s,t)dt. Aix`o ho fem ja que ser`a necessari per a l’argument que emprarem a continuaci´o per a obtenir una equaci´o que ens doni el preu a aquest tipus d’opcions. Com acabem de comentar el valor de la nostra opci´o asi`atica (i m´es generalment d’una que dep`en d’una funci´o f(s,t) sobre el subjacent) ´es una funci´o de tres variables; s,ti la variable d’estat I,V(s,t,I). Considerem ara un argument similar al fet servir anteriorment i ens constru¨ım una cartera amb una opci´o i venem un nombre ∆ d’unitats del subjacent. Π = V(s,t,I)−∆s Aleshores derivant de manera similar a la que hem emprat abans, per`o tenint en compte que ara Vdep`en de Iobtenim: dΠ = ∂V ∂t+1 2σ2s2∂2V ∂2s2dt +∂V ∂IdI +∂V ∂s−∆ds Si prenem ∆ = ∂V ∂sfem una cobertura del risc i resulta dΠ = ∂V ∂t+1 2σ2s2∂2V ∂2s2+f(s,t)∂V ∂Idt Com que sabem que amb la tria ∆ = ∂V ∂sla nostra cartera ´es lliure de risc tenim dΠ = rΠdt i per tant si ho ajuntem tot obtenim l’equaci´o per a modelar el preu de les opcions. ∂V ∂t+1 2σ2s2∂2V ∂2s2+f(s,t)∂V ∂I+rs ∂V ∂s−rV = 0 (19) En el cas d’una opci´o asi`atica amb mitjana aritm`etica cont´ınua obtenim doncs ∂V ∂t+1 2σ2s2∂2V ∂2s2+s∂V ∂I+rs ∂V ∂s−rV = 0 Aquesta equaci´o diferencial estoc`astica tindr`a com a condicions inicials la funci´o de payoff φ(s,I). Per tant, el problema de les opcions asi`atiques es pot resumir en: 1. I(T) = RT 0Sτdτ 2. ∂V ∂t+1 2σ2s2∂2V ∂2s2+s∂V ∂I+rs ∂V ∂s−rV = 0 15 Edps a les finances 3. M`etodes num`erics per l’equaci´o de Black-Scholes 3.1 M`etodes de difer`encies finites La primera fam´ılia de m`etodes num`erics que implementarem es basa en la t`ecnica de les diferencies finites i ens servir`a per a resoldre les opcions europees i tamb´e per a les americanes. Els m`etodes de diferencies finites no s´on tant adients per a resoldre opcions asi`atiques i per a resoldre-les emprarem m`etodes de Montecarlo. Suposem que volem resoldre l’equaci´o de Black-Scholes per a una call amb K= 100 i T= 3 i no ens preocupem de les altres constants. Per tant, considerem resoldre l’equaci´o en la regi´o [0, 3K]×[0, 3]. Ara dividim l’interval [0, 3K] en mintervals de longitud δsi de la mateixa manera partim l’interval [0, 3] en n intervals de longitud δtobtenint una malla de (m+ 1)(n+ 1) nodes. La nostra idea ´es aproximar la funci´o V(s,t) en els nodes de la malla i aproximar els altres punts interpolant amb els punts de la malla. Figura 7: Malla m`etode diferencies finites Numerem els nodes de la seg¨uent forma Vk ion i∈[0, m] i k∈[0, n]. Aleshores tenim Vk i= V(iδs,T−kδt). Comptem el temps a la inversa ja que per k= 0 es dona t=T, que ´es on imposarem la nostra condici´o inicial. 3.1.1 Aproximaci´o de les derivades de V(s,t) Com a l’equaci´o de Black-Scholes apareixen derivades de primer ordre respecte el temps i de primer i segon ordre respecte l’espai hem de construir unes aproximacions num`eriques a partir dels punts de la malla. Prenent la definici´o estricta de la primera derivada respecte el temps tenim ∂V ∂t(s,t) = lim h→0 V(s,t+h)−V(s,t) h i per tant naturalment, usant els punts de la malla, podem dir que ∂V ∂t(s,t)≃Vk i−Vk+1 i δt 16 Fent l’expansi´o de Taylor de V(s,t) podem observar que t´e un error del ordre de O(δt). Per a calcular la primera derivada respecte l’espai usarem diferencies centrades, ´es a dir, en comptes de ∂V ∂s(s,t)≃Vk i+1 −Vk i δs usarem ∂V ∂s(s,t)≃Vk i+1 −Vk i−1 2δs La segona de les dues aproximacions ´es l’obtinguda per el m`etode de les difer`encies centrades i la preferim per el seu menor error. Mentre que la primera de les aproximacions difer`encia endevant t´e un error del ordre de O(δs) la segona t´e un error de O(δs2). Finalment per a poder calcular la segona derivada respecte al temps usem la seg¨uent aproximaci´o ∂2V ∂s2(s,t)≃Vk i+1 −2Vk i+Vk i−1 δs2 que tamb´e t´e error d’ordre O(δs2). 3.2 Opcions europees 3.2.1 M`etode Expl´ıcit El primer dels m`etodes num`erics que implementarem ser`a un model de difer`encies finites expl´ıcit. Si prenem l’equaci´o de Black-Scholes en el tic de temps ktenim ∂V ∂t+1 2σ2s2∂2V ∂2s2+rs ∂V ∂s−rV = 0 (20) i hi substitu¨ım les derivades per les aproximacions usant els valors de la malla obtenim Vk i−Vk+1 i δt+1 2σs2 Vk i+1 −2Vk i+Vk i−1 δs2!+rs Vk i+1 −Vk i−1 2δs!−rV k i=O(δt,δs2) (21) Si reordenem aquesta equaci´o podem a¨ıllar Vk+1 icom a funcions de Vk i−1,Vk iiVk i+1. Vk+1 i=Ak iVk i−1+ (1 + Bk i)Vk i+Ck iVk i+1 (22) on 1. Ak i=δt δs21 2σ2s2−1 2 δt δsrs 2. Bk i=−2δt δs21 2σ2s2−δtr 3. Ck i=δt δs21 2σ2s2+1 2 δt δsrs Observaci´o 3.1.Tant en aquest cas com en el cas del m`etode impl´ıcit i en de Crank-Nicolson no ´es molt important si els coeficients A,B i C s´on avaluats en el tic de temps ko en el tic de temps k+ 1 ja que l’ordre de converg`encia dels m`etodes no es veur`a afectat. D’aquesta manera podem calcular expl´ıcitament el vector Vk+1 (que t´e per components els elements Vk+1 i) a partir de Vk(amb components Vk i). Hi ha un detall, amb aquest m`etode podem calcular tots els elements de Vk+1, excepte per el primer element i tamb´e per a l’´ultim, que els calcularem ajudant-nos de les condicions de contorn. 17 Edps a les finances Condicions de contorn Podem triar entre diverses condicions de contorn. N’expliquem dues de senzilles. Podrem implementar a cada extrem un tipus de soluci´o diferent, no cal usar el mateix m`etode per als dos extrems. 1. Fixem els valors Vk+1 0iVk+1 m. Per un call fixarem Vk+1 0= 0 i Vk+1 m=K−Ke−rδtk . D’altra banda per una put fixarem Vk+1 0=Ke−rδtk iVk+1 m= 0. 2. Podem fixar que la soluci´o tingui derivada segona zero en l’extrem no nul, sempre que l’opci´o tingui un payoff com a molt lineal respecte el subjacent. Usant difer`encies centrades podem obtenir Vk+1 0= 2Vk+1 1−Vk+1 2 i d’aquesta manera un cop coneguts Vk+1 1iVk+1 2a partir del m`etode expl´ıcit ja podem obtenir tots els valors de Vk+1. Aquesta condici´o de contorn t´e un avantatge addicional i es que ´es independent de l’opci´o triada sempre que el contracte tingui un payoff lineal respecte el subjacent. Una vegada implementades les condicions de contorn ja podrem completar el nostre m`etode. L’anomenarem el m`etode expl´ıcit de difer`encies finites per a l’equaci´o de Black-Scholes. Observaci´o 3.2.En aquest esquema a cada pas estem ignorant un error d’ordre O(δt,δs2), anomenant error local de truncament. Aquest m`etode t´e l’avantatge que podem calcular expl´ıcitament la seg¨uent iteraci´o a partir del vector anterior de manera que la seva implementaci´o ´es prou simple. Converg`encia del m`etode i regi´o d’estabilitat Analitzem la converg`encia del m`etode tant en espai com en temps. Ens interessa calcular l’ordre de converg`encia tant en espai com en temps del nostre m`etode. Ho fem calculant el logaritme del error respecte el logaritme del pas de temps o espai en cada cas. Fent els c`alculs obtenim el seg¨uent: calculant el pendent de les dues rectes podem observar la converg`encia d’ordre (a) Converg`encia en l’espai (b) Converg`encia en el temps Figura 8: Gr`afiques de converg`encia per el m`etode expl´ıcit 2 en l’espai i d’ordre 1 en el temps. El major desavantatge del m`etode expl´ıcit ´es la seva inestabilitat, ´es a dir, es pot provar que si δt≤δs2 σ2s2 18 no es compleix aleshores el m`etode no ser`a estable i la seva converg`encia no estar`a garantida. Aix`o ens indica (prenent el denominador com a constant) que per a reduir el pas d’espai a la meitat hem de fer el pas de temps quatre vegades m´es petit, la limitaci´o m´es important del m`etode. 3.2.2 M`etode Impl´ıcit Usem de la mateixa manera que amb el m`etode expl´ıcit l’equaci´o de Black-Scholes, per`o ara en el tic de temps k+ 1. Obtenim doncs la seg¨uent equaci´o: Vk i−Vk+1 i δt+1 2σs2 Vk+1 i+1 −2Vk+1 i+Vk+1 i−1 δs2!+rs Vk+1 i+1 −Vk+1 i−1 2δs!−rV k+1 i=O(δt,δs2) (23) Ara podem a¨ıllar Vk ii procedint de la mateixa manera que en el cas expl´ıcit per`o situats en el tic de temps k+ 1. Obtenim doncs el seg¨uent esquema amb el que podem calcular Vk ia partir de Vk+1 i−1,Vk+1 iiVk+1 i+1 . Vk i=Ak+1 iVk+1 i−1+ (1 + Bk+1 i)Vk+1 i+Ck+1 iVk+1 i+1 (24) on 1. Ak i=−δt δs21 2σ2s2−1 2 δt δsrs 2. Bk i= 2 δt δs21 2σ2s2+δtr 3. Ck i=−δt δs21 2σ2s2+1 2 δt δsrs Usant aquestes equacions podem construir un sistema d’equacions la soluci´o del qual sigui el vector Vk+1.          Ak+1 21 + Bk+1 2Ck+1 20··· 0 0Ak+1 31 + Bk+1 3Ck+1 3 ...0 . . .............. . . . . ....Ak+1 m−21 + Bk+1 m−2Ck+1 m−20 0··· 0Ak+1 m−11 + Bk+1 m−1Ck+1 m−1                     Vk+1 1 Vk+1 2. . . . . . . . . Vk+1 m            =            Vk 2 Vk 3 . . . . . . . . . Vk m−1            (25) Aix`o ens indica que podem calcular l’iteraci´o k+ 1 a partir de Vki una matriu que ens acabar`a resultant constant en el temps en el nostre cas. La nostra matriu tridiagonal tindr`a un tamany (m−2)mde manera que calen dues condicions m´es per a tal de poder resoldre el sistema de manera ´unica. Aquest ´es el desavantatge del m`etode impl´ıcit, per cada iteraci´o hem de resoldre un sistema d’equacions lineal i no podem calcular-la de manera expl´ıcita, fent que requereixi m´es c`alculs per iteraci´o. Condicions de contorn Hem de tenir en compte ara les condicions de vora. En podem considerar de diferents tipus de la mateixa manera que en el m`etode expl´ıcit, per`o amb alguns petits detalls. 19 Edps a les finances 1. Si volem fixar els valors Vk+1 0iVk+1 mreescrivim el nostre sistema d’equacions lineals eliminant la primera i l’´ultima columna (les que afecten sobre els punts de la vora) de la nostra matriu i afegint els termes que s’eliminen com a un vector per tal de fer que el nostre sistema tingui soluci´o ´unica al tenir la nostra matriu mida (m−2)(m−2).          Bk+1 2Ck+1 20··· 0 Ak+1 31 + Bk+1 3Ck+1 3 ...0 . . ........... . . . . .Ak+1 m−21 + Bk+1 m−2Ck+1 m−20 0··· Ak+1 m−11 + Bk+1 m−1Ck+1 m−1                     Vk 2 Vk 3 . . . . . . . . . Vk m−1            =            Vk 2 Vk 3 . . . . . . . . . Vk m−1            −           Ak+1 2Vk+1 1 0 . . . . . . 0 Ck+1 m−1Vk+1 m           2. Si volem fixar la segona derivada com a nul·la hem d’incorporar la condici´o Vk+1 0= 2Vk+1 1−Vk+1 2 al nostre sistema d’equacions. Ho fem substituint Mk+1Vk+1 per          1 + Bk+1 2+ 2Ak+1 2Ck+1 2−Ak+1 20··· Ak+1 31 + Bk+1 3Ck+1 3 ... ............ ...Ak+1 m−21 + Bk+1 m−2Ck+1 m−2 ··· 0Ak+1 m−11 + Bk+1 m−1                     Vk 1 Vk 2 . . . . . . . . . Vk m            on Mk+1 ´es la matriu del nostre sistema en (25). Observaci´o 3.3.Podem usar en cada un dels extrems cada unes de les t`ecniques descrites, no cal usar per els dos extrems el mateix tipus de condicions de vora. Converg`encia del m`etode i regi´o d’estabilitat Analitzem ara l’ordre de converg`encia del nostre m`etode la mateixa manera que ho vam fer per a la variant expl´ıcita. Obtenim els mateixos ordres de converg`encia (a) Converg`encia en l’espai (b) Converg`encia en el temps Figura 9: Gr`afiques de converg`encia per el m`etode impl´ıcit que en el m`etode expl´ıcit. S’observa una converg`encia d’ordre 2 respecte l’espai i d’ordre 1 respecte el 20 temps. El principal avantatge respecte el m`etode expl´ıcit recau en que ´es un m`etode que sempre ´es num`ericament estable, sense importar la relaci´o entre δtiδs. Aix`o ens permetr`a evitar haver de triar un δtmolt petit si volem fer prou petit δs. 3.2.3 M`etode Crank-Nicolson El m`etode de Crank-Nicolson ´es una combinaci´o dels dos m`etodes que acabem de descriure. B`asicament per a cada iteraci´o consisteix en calcular la soluci´o amb la mitjana aritm`etica entre el m`etode expl´ıcit i el m`etode impl´ıcit. Vk i−Vk+1 i δt+1 4σs2 Vk+1 i+1 −2Vk+1 i+Vk+1 i−1 δs2!+1 4σs2 Vk i+1 −2Vk i+Vk i−1 δs2! +rs 2 Vk+1 i+1 −Vk+1 i−1 2δs!+rs 2 Vk i+1 −Vk i−1 2δs!−rV 2Vk+1 i−rV 2Vk i=O(δt,δs2) (26) En aquesta primera equaci´o estem fent la mitjana aritm`etica de les equacions (20)i(23) corresponents al m`etode expl´ıcit i impl´ıcit. En aquest cas no a¨ıllarem Vk isin´o que amb construir un sistema de la mateixa forma que en els dos m`etodes anteriors (expressant l’equaci´o com a suma dels nodes implicats multiplicats per els seus factors corresponents) n’hi haur`a suficient. Obtenim aquest esquema −Ak+1 iVk+1 i−1+ (1 −Bk+1 i)Vk+1 i−Ck+1 iVk+1 i+1 =Ak iVk i−1+ (1 + Bk i)Vk i+Ck iVk i+1 amb aquests coeficients: 1. Ak i=1 2 δt δs21 2σ2s2+1 4 δt δsrs 2. Bk i=−δt δs21 2σ2s2+1 2δtr 3. Ck i=1 2 δt δs21 2σ2s2−1 4 δt δsrs D’aquestes equacions obtenim un sistema d’equacions lineal de la forma Mk+1Vk+1 +rk=NkVk.          Ak+1 21 + Bk+1 2Ck+1 20··· 0 0Ak+1 31 + Bk+1 3Ck+1 3 ...0 . . .............. . . . . ....Ak+1 m−21 + Bk+1 m−2Ck+1 m−20 0··· 0Ak+1 m−11 + Bk+1 m−1Ck+1 m−1                     Vk+1 1 Vk+1 2. . . . . . . . . Vk+1 m            =          −Ak 21−Bk 2−Ck 20··· 0 0−Ak 31−Bk 3−Ck 3 ...0 . . .............. . . . . ....−Ak m−21−Bk m−2−Ck m−20 0··· 0−Ak m−11−Bk m−1−Ck m−1                     Vk 1 Vk 2 . . . . . . . . . Vk m            (27) 21 Edps a les finances Condicions de contorn De manera similar a com es va fer amb el m`etode impl´ıcit hem d’implementar les condicions de vora per aconseguir que el sistema tingui soluci´o ´unica. Tenim diverses opcions: •Si volem fixar els valors de Vk+1 1iVk+1 maleshores usem la mateixa t`ecnica que en el cas del m`etode impl´ıcit sobre la matriu Mk+1. N’eliminem la primera i la ´ultima de les columnes i afegim els termes eliminats en forma de vector lliure. Matricialment substitu¨ım Mk+1Vk+1 per          1 + Bk+1 2Ck+1 20··· Ak+1 31 + Bk+1 3Ck+1 3 ... ............ ...Ak+1 m−21 + Bk+1 m−2Ck+1 m−2 ··· 0Ak+1 m−11 + Bk+1 m−1                     Vk 1 Vk 2 . . . . . . . . . Vk m            +           Ak+1 2Vk+1 1 0 . . . . . . 0 Ck+1 m−1Vk+1 m           •Si volem fixar la segona derivada com a nul·la hem d’incorporar la condici´o Vk+1 0= 2Vk+1 1−Vk+1 2 al nostre sistema d’equacions. Ho fem de la mateixa manera que en m`etode impl´ıcit, substituint Mk+1Vk+1 per          1 + Bk+1 2+ 2Ak+1 2Ck+1 2−Ak+1 20··· Ak+1 31 + Bk+1 3Ck+1 3 ... ............ ...Ak+1 m−21 + Bk+1 m−2Ck+1 m−2 ··· 0Ak+1 m−11 + Bk+1 m−1                     Vk 1 Vk 2 . . . . . . . . . Vk m            Converg`encia del m`etode i regi´o d’estabilitat Si estudiem l’ordre de converg`encia del nostre m`etode respecte l’espai i del temps de la mateixa manera que l’explicada en el m`etode expl´ıcit podem obtenir una converg`encia d’ordre 2 tant en temps com en espai. Aix`o es una millora respecte els m`etodes expl´ıcit i impl´ıcit que nom´es tenen ordre 1 de converg`encia respecte el temps. Aquesta ´es un dels grans avantatges del m`etode de Crank-Nicolson i el perqu`e a la pr`actica quasi mai es fa servir el m`etode impl´ıcit ja que si per a cada iteraci´o hem de resoldre un sistema d’equacions lineals millor fer servir el m`etode de Crank-Nicolson que ens assegura un millor ordre de converg`encia. Respecte a la regi´o d’estabilitat del m`etode es pot demostrar que el m`etode ´es sempre estable, sense importar la tria de δtiδs. 3.2.4 Resoluci´o de sistemes lineals Tant en el m`etode impl´ıcit com en el m`etode de Crank-Nicolson necessitem resoldre sistemes lineals a cada pas per a trobar Vk+1. Normalment ho fem amb la rutina que ja t´e Matlab, per`o per a l’implementaci´o a les opcions americanes necessitarem tenir implementat el m`etode de sobre-relaxaci´o successiva. Suposem que volem resoldre un sistema est`andard d’equacions lineal Ax =bon A ´es una matriu n×n d’elements aij i tant xcom bs´on vectors de mida n. Primerament partim Aen suma de tres matrius A= D+L+Uon D´es la seva diagonal i L,Ules seves matrius triangulars inferiors i superiors respectivament. D’aquesta manera podem reescriure el sistema d’equacions d’aquesta manera: (D+ωL)x=ωb−(ωU+ (ω−1)D)x 22 Com podem observar en la Figura 14 es pot observar que el preu d’una opci´o asi`atica ´es una mica diferent Figura 14: En aquesta figura veiem representats els preus de tres tipus d’opcions sota els mateixos par`ametres r= 0.08, σ= 0.25, K= 100 i m= 100 per a l’opci´o asi`atica. del de la seva contrapartida europea. Aix`o ´es aix´ı ja que la opci´o europea nom´es t´e en compte el preu a temps final mentre que l’asi`atica t´e en compte tota l’evoluci´o del preu, cosa que la protegeix una mica m´es contra la volatilitat del preu. Tamb´e podem comparar l’efecte dels par`ametres que podem observar en la figura respecte riσcom podem veure en la Figura 15. Com podem veure l’augment de la volatilitat fa augmentar el preu de les opcions, en la figura nom´es apareixen opcions de venda per`o augmenta el preu de totes les opcions. D’altra banda un augment en el tipus d’inter`es augmenta el valor de les opcions de compra ja que fa augmentar el preu final esperat, i per aix`o fa disminuir de valor les opcions de venda. 29 Edps a les finances (a) Preu d’una opci´o europea per diferents valors de r(b) Preu d’una opci´o europea per diferents valors de σ Figura 15: En aquestes dues figures representem els preus d’opcions de compra europees respectivament per als par`ametres seg¨uents: r= 0.08, T= 3, K= 100,i per a diversos valors de σ. Apareixen tamb´e en blau les liquidacions de les opcions. 5. Bibliografia Refer`encies [1] Karel in ’t Hout. Numerical Partial Differential Equations in Finance Explained, Palgrave Macmillan, 2017. [2] Paul Wilmott. Derivatives. The theory and practice of financial engineering, John Wiley & Sons, 1998. [3] Paul Wilmott, Jeff Dewynne, Sam Howison. Option pricing. Mathematical models and computation, Oxford Financial Press, 1994. [4] Hongbin Zhang. Pricing Asian Options using Monte Carlo Methods, Project Report Uppsala Universitet, 2009. 30 A. Codis M`etode Diferencies Finites A.1 M`etode Expl´ıcit function U = p a r a b o l i c E u l e r ( x , Ax , At , nOfSteps , f , sigma , r , option ,K) %x : d i s c r e t i s a t i o n of the S i n t e r v a l % Ax At space / time s t e p s % nOfSteps : number of time s t e p s % f : i n i t a l c o n d i t i o n ( p a yo ff ) % sigma r K, o ption parameters % option = 0 c a l l // option = 1 put NAssetSteps = size( f ) ; %Initialization Un = zeros (NAssetSteps ); Derivative = zeros ( NAssetSteps ) ; SDerivative = zeros (NAssetSteps ); Unm1 = f ; %s o l u t i o n at time n−1 U = f ; ha lf s qs ig ma = (1/2)∗sigma ˆ2; xsqd = x . ∗x ; %Loop i n time s t e p s for n=1:nOfSteps for i = 2:NAssetSteps(2)−1 %We approximate d e r i v a t i v e and 2n d e r i v a t i v e D e r i v a t i v e ( i ) = (Unm1( i + 1) −Unm1( i −1))/(2∗Ax ) ; S D e r i v a t i v e ( i ) = (Unm1( i + 1) −2∗Unm1( i ) + Unm1( i −1)) /( Ax ˆ 2) ; %i n t e r i o r node Un( i ) = Unm1( i ) + At ∗(halfsqsigma∗xsqd ( i )∗SDerivative(i) + r∗x ( i )∗Derivative(i) −r∗Unm1( i ) ) ; end %Boundary c o n d i t i o n s i f option == 0 Un(1) = 0 ; Un( end ) = 2∗Un(end −1) −Un( end −2 ) ; else Un(1) = K∗exp(−r∗At∗n ) ; Un( end ) = 0 ; end Unm1 = Un ; %next step U = Un ; end 31 Edps a les finances A.2 M`etode Impl´ıcit function U = i m p l i c i t E u l e r (x , Ax , At , nOfSteps , f , sigma , r , option ,K) %x : d i s c r e t i s a t i o n of the S i n t e r v a l % Ax At space / time s t e p s % nOfSteps : number of time s t e p s % f : i n i t a l c o n d i t i o n ( pa y of f ) % sigma r K, o ption parameters % option = 0 c a l l // option = 1 put ha lf s qs ig ma = (1/2)∗sigma ˆ2; function r e s = A( n ) r e s = ( (1/2 )∗r .∗n−halfsqsigma.∗n .∗n)∗At ; end function r e s = B( n ) r e s = ( sigma ˆ2.∗n .∗n + r )∗At ; end function r e s = C( n ) r e s = −(halfsqsigma .∗n .∗n + (1/2)∗r .∗n )∗At ; end NAssetSteps = size( f ) ; %Initialization Unm1 = f ; %s o l u t i o n at time n−1 U = f ; %s o l u t i o n v e ct o r ( i n i t i a t e d by IC ) %Matrix Assembly % d i a g o n a l s MA = A( 0 : NAssetSteps (2) −1 ) ; M A(1) = [ ] ; M A(1) = [ ] ; M B = B( 0 : NAssetSteps (2) −1) + ones ( NAssetSteps ) ; M B(1) = [ ] ; M B( end ) = [ ] ; M C = C( 0 : NAssetSteps (2) −1 ) ; M C( end ) = [ ] ; M C( end ) = [ ] ; %s p a r s e matr ix c r e a t i o n M = spdiags ( [ M A;M B;M C] ’ , −1:1,NAssetSteps(2) −2,NAssetSteps(2) −2 ) ; %matrix d i s p l a y %[M A;M B; M C] %f u l l (M) %Loop i n time s t e p s for n=1:nOfSteps F = Unm1; F(1) = [ ] ; F( end ) = [ ] ; 32 i f option == 0 %boundary c o n d i t i o n s U 0 = 0; F(1) = F(1) −A(1)∗U 0 ; U end = x ( end )−K∗exp(−r∗At ∗(n ) ) ; F( end ) = F(end )−C(NAssetSteps(2) −2)∗U end ; %s o l u t i o n of the l i n e a l system U = M\F ’ ; end i f option == 1 %boundary c o n d i t i o n s U 0 = K∗exp(−r∗At ∗( n ) ) ; F(1) = F(1) −A(1)∗K∗exp(−r∗At ∗( n ) ) ; F( end ) = F(end )−C(NAssetSteps (2))∗0 ; U end = 0; %s o l u t i o n of the l i n e a l system U = M\F ’ ; end U = [ U 0 U’ U end ] ; Unm1 = U; %next step end end A.3 M`etode Crank-Nicolson f u n c t i o n U = crank ( x , Ax , At , nOfSteps , f , sigma , r , option ,K) ha lf s qs ig ma = (1/2)∗sigma ˆ2; %Wilmott pg638 f u n c t i o n s f u n c t i o n r e s = A( n ) r e s = (1/ 2)∗((1/2)∗r .∗n−halfsqsigma.∗n .∗n)∗At ; end f u n c t i o n r e s = B( n ) r e s = ( 1/2) ∗( sigma ˆ2.∗n . ∗n + r )∗At ; end f u n c t i o n r e s = C( n ) r e s = ( h al fs qs i gm a . ∗n .∗n + (1/2)∗r .∗n )∗At ∗(−1/2); end NAssetSteps = s i z e ( f ) ; %Initialization Unm1 = f ; %s o l u t i o n at time n−1 U = f ; %s o l u t i o n v ec t or ( i n i t i a t e d by IC ) 33 Edps a les finances xsqd = x . ∗x ; %Matrix Assembly % d i a g o n a l s M A = A( 0 : NAssetSteps (2) −1 ) ; M A(1) = [ ] ; M A(1) = [ ] ; M B = B( 0 : NAssetSteps (2) −1) + ones ( NAssetSteps ) ; M B(1) = [ ] ; M B( end ) = [ ] ; M C = C( 0 : NAssetSteps (2) −1 ) ; M C( end ) = [ ] ; M C( end ) = [ ] ; N A = −A(0: NAssetSteps(2) −1 ) ; N A(1) = [ ] ; N A( end ) = [ ] ; N B = −B(0: NAssetSteps(2) −1) + ones ( NAssetSteps ) ; N B(1) = [ ] ; N B( end ) = [ ] ; N C = −C(0: NAssetSteps (2) −1 ) ; N C (1) = [ ] ; N C( end ) = [ ] ; %s p a r s e matr ix c r e a t i o n M = s pd ia g s ( [ M A;M B;M C] ’ , −1:1,NAssetSteps(2) −2,NAssetSteps(2) −2 ) ; N = s pd ia gs ( [ N A ; N B ; N C ] ’ , 0 : 2 , NAssetSteps (2)−2,NAssetSteps (2)); %boundary c o n d i t i o n s M(1 ,1) = M(1 ,1) + 2∗A( 1 ) ; M(1 ,2) = M(1 ,2) −A( 1 ) ; M( end , end ) = M( end , end ) + 2∗C(NAssetSteps(2) −2 ) ; M( end , end −1) = M( end , end −1) −C(NAssetSteps(2) −2 ) ; %matrix d i s p l a y %[M A;M B;M C] %f u l l (M) %f u l l (N) %Loop i n time s t e p s f o r n=1: nOfSteps F = Unm1; F = N∗F ’ ; %s o l u t i o n of the l i n e a l system U = M\F ; U 0 = 2∗U(1) −U(2); U end = 2∗U( end ) −U( end −1 ) ; U = [ U 0 U’ U end ] ; Unm1 = U; %next s tep end end 34 A.4 M`etode opcions americanes A.5 M`etode de sobre-relaxaci´o successiva per a sistemes lineals f u n c t i o n [U, omega , optimum reached , n I t s ] = SOR( Lo ,D, Up , q , n , option ,K, x , omega , optimum reached , n I t s ) %SOR s o l u t i o n ( wilmott 646) t o l = 1e −14; v = z e r o s ( 1 , n ) ; e r r o r = 1; I t s = 0; wh ile e r r o r >tol e r r o r = 0; f o r i = 2: n −1 i f opt ion == 0 temp = max( v ( i ) + ( omega/D( i ) ) ∗ (q ( i −1) −Up( i ) ∗v ( i + 1) −D( i ) ∗v ( i ) −Lo ( i ) ∗v ( i −1) ) , x ( i ) −K) ; e l s e temp = max( v ( i ) + ( omega/D( i ) ) ∗ (q ( i −1) −Up( i ) ∗v ( i + 1) −D( i ) ∗v ( i ) −Lo ( i ) ∗v ( i −1) ) , K−x ( i ) ) ; end e r r o r = e r r o r + ( temp −v ( i ) ) ∗( temp −v ( i ) ) ; v ( i ) = temp ; end I t s = I t s + 1 ; end v (1) = [ ] ; v ( end ) = [ ] ; U = v ’ ; if optimum reached == false i f I t s <n I t s omega = omega + 0 . 0 5 ; n I t s = I t s ; e l s e optimum reached = true ; end end end 35 Edps a les finances B. Codi M`etode de Montecarlo B.1 M`etode opcions europees n = 10000; S 0 = 100; K = 100; r = 0 . 0 8 ; sigma = 0 . 2 5 ; T = 3; opti on = 0; %0 = c a l l / 1 = put %ANTITHETIC VARIATE METHOD sum = 0; f o r i = 0: n rand = normrnd ( 0 , 1 ) ; S T = S 0 ∗exp (( r −( sigma ˆ2/2))∗T + sigma ∗s q r t (T)∗rand ) ; S T2 = S 0∗exp ( ( r −( sigma ˆ2/2))∗T−sigma ∗s q r t (T)∗rand ) ; i f opt ion == 0 sum = sum + max( S T −K, 0 ) + max( S T2 −K,0); e l s e sum = sum + max(K −S T , 0 ) + max(K −S T2 , 0 ) ; end end r e s = exp(−r∗T)∗( sum/(2∗n ) ) %CRUDE MONTECARLO METHOD sum = 0; f o r i = 0: n S T = S 0 ∗exp (( r −( sigma ˆ2/2))∗T + sigma ∗s q r t (T)∗normrnd ( 0 , 1 ) ) ; i f opt ion == 0 sum = sum + max( S T −K, 0 ) ; e l s e sum = sum + max(K −S T , 0 ) ; end end r e s = exp(−r∗T)∗( sum/n ) B.2 M`etode opcions asi`atiques n = 10000; m = 100; S 0 = 100; 36 K = 100; r = 0 . 0 8 ; sigma = 0 . 2 5 ; T = 3; opti on = 0; %0 = c a l l / 1 = put %ANTITHETIC VARIATE METHOD sum = 0; f o r i = 0: n S T = S 0 ; S T2 = S 0 ; parsum1 = 0; parsum2 = 0; f o r j = 0:m rand = normrnd ( 0 , 1 ) ; S T = S T∗exp (( r −( sigma ˆ2/2))∗(T/m) + sigma ∗s q r t (T/m)∗rand ) ; S T2 = S T2∗exp (( r −( sigma ˆ2/2))∗(T/m) −sigma ∗s q r t (T/m)∗rand ) ; parsum1 = parsum1 + S T ; parsum2 = parsum2 + S T2 ; end i f opt ion == 0 sum = sum + max( parsum1 /(m + 1) −K, 0 ) + max( parsum2 /(m + 1) −K , 0 ) ; e l s e sum = sum + max(K −parsum1 /(m + 1) ,0) + max(K −parsum2 /(m + 1 ) , 0 ) ; end end r e s = exp(−r∗T)∗( sum/(2∗n ) ) %CRUDE MONTECARLO METHOD sum = 0; f o r i = 0: n ST = S 0 ; parsum = 0; f o r j = 0:m S T = S T∗exp (( r −( sigma ˆ2/2))∗(T/m) + sigma ∗s q r t (T/m)∗normrnd ( 0 , 1 ) ) ; parsum = parsum + S T; end i f opt ion == 0 sum = sum + max( parsum /(m + 1) −K, 0 ) ; e l s e sum = sum + max(K −parsum /(m + 1 ) , 0 ) ; end end r e s = exp(−r∗T)∗( sum/n ) 37