Universit Paris I Panthon-Sorbonne

Master 2 MAEF - Spcialit Recherche : MMMEF

Ecole Nationale Suprieure de Techniques Avances

Ismail LAACHIR
Mthode des Points Intrieurs en Optimisation et application la calibration du modle d'Heston

&

Rsum : La pertinence d'un modle est souvent value par plusieurs critres parmi lesquels la capacit des prix des options vanilles gnrs par le modle coller aux prix du march. En eet, le but de la calibration est de pouvoir utiliser le modle, une fois calibr, pour calculer le prix des options exotiques et la stratgie de leur couverture avec ces options vanilles. Cela impose de pouvoir choisir les paramtres du modle de manire retrouver les volatilits implicites de march. Cette partie essentielle correspond la calibration. Il est donc indispensable d'utiliser un algorithme d'optimisation qui soit rapide, prcis et surtout numriquement stable. Dans ce sens, la premire partie de mon stage a consist en l'implmentation d'un algorithme d'optimisation bas sur la mthode des Points Intrieurs propos par Armand et al [1]. Cet algorithme a ensuite t test sur plusieurs problmes d'optimisation classiques pour valuer la convergence et la robustesse de l'algorithme. Ensuite, dans la deuxime et troisime parties, nous avons utilis cet algorithme pour la calibration d'un modle volatilit stochastique qu'est le modle d'Heston, aprs une tude pralable du modle et du rle de chacun de ses paramtres.

Rapport de stage prsent le 15 septembre 2008.

Jury :
David LEFEVRE

Enseignant-chercheur ENSTA UMA/OC (Pice 3055) 32, boulevard Victor F-75379 Paris cedex 15
Georges NEMES

Socit Gnrale Direction des Risques RISQ/RDM/MOD

Table des matires

Introduction 1 Mthode des Points Intrieurs en Optimisation : Algorithme BFGS

2 4

1.1 1.2

1.3

Introduction . . . . . . . . . . . . . . . . . . . . Algorithme BFGS d'Armand et al [1] . . . . . . 1.2.1 Calcul d'une direction de descente . . . 1.2.2 Recherche Linaire . . . . . . . . . . . . 1.2.3 Actualisation de l'information du second 1.2.4 Algorithmes A et A . . . . . . . . . . . 1.2.5 Convergence de l'algorithme global A . Rsultats Numriques . . . . . . . . . . . . . . 1.3.1 Fonction de Rosenbrock . . . . . . . . . 1.3.2 Welded Beam Design Problem . . . . . 1.3.3 Pressure Vessel Design Problem . . . . .

. . . . . . . . . . . . . . . . ordre . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 4 6 6 7 8 10 10 11 11 12
15

Modle d'Heston : Evaluation

2.1 2.2

Prsentation du modle . . . . . . . . . . . . . . . . 2.1.1 Proprits thoriques . . . . . . . . . . . . . Prix d'un Call : Formule semi-explicite . . . . . . . . 2.2.1 Instabilit numriques : Logarithme complexe 2.2.2 Intgration Numrique . . . . . . . . . . . . .

15 15 16 17 19
20

Modle d'Heston : Calibration

3.1 3.2

Volatilit Implicite : Eet des paramtres . . . . . . . . . . . . . . Calibration : Introduction . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Application aux prix gnrs par le modle d'Heston . . . . 3.2.2 Calibration du modle d'Heston sur des donnes de march Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A Modle d'Heston : Rsultats thoriques

20 23 23 24 27
28

A.1 Existence du processus (Vt )t0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Contrle uniforme sur les moments de X et V . . . . . . . . . . . . . . . . . . . . . . . A.3 Continuit de f1 et f2 en 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B Problmes test

28 29 30
34

Introduction

Pour qu'un modle puisse tre pertinent dans la pratique, il doit tre capable de redonner les prix de march des options vanilles. Cela impose de pouvoir choisir les paramtres du modle de manire retrouver les volatilits implicites de march. Cette partie essentielle correspond la calibration. La calibration est ce que l'on appelle un problme mal pos (ill-posed problem) : il est quasiment impossible qu'un modle paramtrique soit en mesure de retrouver exactement tous les prix de march. On tente alors de minimiser une fonction erreur qui mesure la proximit de la surface des prix thoriques de celle du march. A l'inverse, cette minimisation peut gnralement avoir plusieurs solutions et dpendre fortement de la fonction erreur que l'on choisit. C'est pour ces raisons qu'il est indispensable d'utiliser un algorithme d'optimisation qui soit rapide, prcis et surtout numriquement stable, dans le sens o un petit changement des prix observs ne doit pas entraner de grands changements des paramtres. Il existe plusieurs types d'algorithmes d'optimisation, parmi lesquels on distingue les mthodes de Points Intrieurs, qui ont suscit un intrt croissant ces deux dernires dcennies. Dans ce sens, la premire partie de mon stage a consist en l'implmentation d'un algorithme d'optimisation bas sur la mthode des Points Intrieurs propos par Armand et al [1]. Cet algorithme a ensuite t test sur plusieurs problmes d'optimisation classiques pour valuer la convergence et la robustesse de l'algorithme. Ces problmes test montrent en eet que l'algorithme est assez prometteur pour son utilisation en calibration. La deuxime partie de ce rapport est consacre un modle volatilit stochastique qu'est le modle d'Heston. Le choix de ce modle est motiv par les constatations exprimentales sur les courbes de volatilit implicite remarques sur le march et les limites des modles jusque l utiliss pour capturer cette dynamique du smile. Ce modle utilise plusieurs paramtres ce qui en fait un modle exible et pouvant gnrer plusieurs formes de smile remarques sur la march. Ensuite, nous avons rappel la mthode d'valuation des options vanilles dans ce modle et les problmes numriques rencontrs lors de son utilisation et les solutions apportes ces problmes pour pouvoir valuer les prix de faon stable et robuste numriquement. Ceci est indispensable pour l'tape de la calibration, qui est tudie dans la troisime partie du rapport. Pour cela nous avons commenc par prsenter la problmatique de la calibration, puis tudi l'eet des dirents paramtres du modle d'Heston sur les volatilits implicites. Ensuite nous avons soumis l'algorithme d'optimisation un dernier test qui est de rcuprer les paramtres avec lesquels une nappe de prix a t gnre. A l'issue des rsultats concluants de ce test, nous avons utilis l'algorithme pour la calibration du modle d'Heston sur des donnes relles de march. La comparaison de nos rsultats avec ceux de l'article [7] montrent la prcision et la robustesse de cet algorithme.

Chapitre 1
Mthode des Points Intrieurs en Optimisation : Algorithme BFGS

1.1 Introduction
L'optimisation est un domaine assez rcent qui a connu un grand essor et gagn en importance dans direntes problmatiques d'ingnierie, de science et d'conomie, protant du dveloppement acclr de la puissance des ordinateurs. Ces dirents besoins ont donn naissance une multitude de branches dans l'optimisation, chacune avec sa thorie et ses mthodes numriques. On peut distinguer par exemple plusieurs dualits : optimisation continue/combinatoire, direntielle/non direntielle, dimension innie/dimension nie et nalement linaire/non linaire. Selon [2] : En optimisation non linaire, on peut ainsi distinguer plusieurs vagues : mthodes de pnalisation, mthode du lagrangien augmente (1958), mthodes de quasi-Newton (1959), mthodes newtoniennes ou SQP (1976), algorithmes de points intrieurs (1984). Une vague n'eace pas la prcdente mais permet d'apporter de meilleures rponses certaines classes de problmes, comme ce fut le cas pour les mthodes de points intrieurs en optimisation semi-dnie positive (SDP). Comme notre besoin est d'avoir un algorithme pour extraire un nombre ni de paramtres, en manipulant des fonctions assez rgulires, nous allons nous intresser un algorithme d'optimisation direntielle non linaire et en dimension nie. Pour notre tude nous allons implmenter un algorithme de Points Intrieurs, et en prenant en compte aussi l'impratif de rapidit d'excution, nous avons choisi d'utiliser un algorithme qui utilise la mthode de quasi-Newton pour calculer l'information de second ordre un moindre eort. Ceci diminuera le nombre d'valuations de la fonction minimiser, ce qui peut se rvler crucial pour l'utilisation pratique de la calibration.

1.2 Algorithme BFGS d'Armand et al [1]

Un problme d'optimisation classique s'crit de la faon suivante :

min f (x) c(x) 0

Hypothse H

(1.1)

Nous supposerons que les fonctions f : Rn R et c : Rn Rm sont direntiables deux fois, f convexe, et les ci concave et qu'au moins une des fonctions f, c1 , ..., cm est fortement convexe 1 . On suppose aussi qu'il existe x tel que c(x) > 0.
1 Une
fonction

: Rn R

est fortement convexe de module

> 0,

si la fonction

(.)

||.|| 2

est convexe.

Ce problme est souvent transform, sous certaines conditions, en un systme d'quations et inquations connu sous le nom des conditions KKT (Karush-Kuhn-Tucker) : (x, ) = f (x) c(x) = 0 C(x) = 0 (P ) c(x) 0 0 o C(x) = diag(c1 (x), ..., cm (x)) et (x, ) = f (x) c(x) reprsente le lagrangien du problme (1.1). La variable Rm est appele vecteur dual, et le systme (P ) primal-dual. La deuxime quation porte le nom des conditions de complmentarit : Elles signient que soit i = 0 soit ci (x) = 0, ce qui reprsente 2m cas tudier. Elle est donc une source principale de dicult pour rsoudre le problme (P). Cette remarque est la base de la mthode des Points Intrieurs. En eet, pour liminer le problme combinatoire de la 2me quation de (P), on peut approximer les conditions d'optimalit exprimes dans (P) par les quations perturbes suivantes : f (x) c(x) = 0 C(x) = e (P ) (c(x), ) > 0 , o est un paramtre strictement positif que l'on fera tendre vers zro. La limite de ces solutions quand tends vers 0 concidera avec la solution de (1.1).
Remarque :

Une deuxime faon, quivalente, d'introduire la mthode des Points Intrieurs est de pnaliser le problme d'optimisation (1.1) comme suit :
min (x) c(x) > 0
m

, avec (x) = f (x)

k=1

log(ci (x)).

Cette fonction, appele "barrire logarithmique", est dnie seulement pour c(x) > 0, ie l'intrieur de l'ensemble des points admissibles, d'o le nom de Points Intrieurs. Le systme d'quations P est alors une formulation des conditions d'optimalit de ce problme pnalis.
L'algorithme BFGS propos dans l'article [1] que nous allons implmenter et utiliser en calibration appartient la classe de mthodes de points intrieurs. Il combine les techniques des Points Intrieurs, et les mthodes de Quasi-Newton pour rsoudre une suite de problmes (P )0 . Cette suite de problmes rsoudre est une dicult supporter pour pouvoir se passer des quations de complmentarit, source principale de complexit. On appelle A l'algorithme pour rsoudre (P ) et A l'algorithme global pour rsoudre P . L'algorithme A peut tre rsum en trois tapes : 1. Calcul d'une direction de descente ie un vecteur d = (dx , d ). 2. Recherche Linaire : chercher un coecient pour que le nouvel itr x+ = x + dx et + = + d vrie un certain critre de descente. 3. Actualisation de l'information du second ordre.

1.2.1 Calcul d'une direction de descente

Les deux quations de (P ) seront rsolues par mthode de quasi-Newton. Cette mthode approxime l'quation vectorielle F (y) = 0, o F : Rp Rp , par une suite de points qui satisfont F (yk ) + F (yk )(yk+1 yk ) = 0, ce qui revient approximer la fonction F par sa tangente en yk . La suite yk est alors dnie par yk+1 = yk F (yk )1 F (yk ), condition que F (yk ) soit inversible. Le vecteur F (yk )1 F (yk ) est appel direction de descente de Newton. Cette mthode est assez simple programmer et assez rapide, surtout quand le gradient de F est facile calculer. Son utilisation pour le problme P donne :

M c(x) c(x) C(x)

, avec = diag(1 , ..., m ) et M = Ce qui donne, en liminant d :
2 xx

dx d (x, ).

f (x) + c(x) e C(x)

(M +

c(x)C(x)1 c(x)T )dx = f (x) + c(x)C(x)1 e = (x)

Il est donc vident que si M est strictement positive, et que (x, ) est strictement admissible ie : c(x) > 0 et > 0, alors la vecteur dx existe et est une direction de descente pour la fonction , ie : (x)dx < 0. Cela nous assure que la fonction (.) diminue au voisinage de x. Par contre, le seul choix xk+1 = xk + dx et k+1 = k + d ne peut garantir lui seul de diminuer la fonction (.). Pour ce but, on utilise souvent la mthode Recherche Linaire.

1.2.2 Recherche Linaire

Pour forcer la diminution de le long de la direction dx , il est possible de calculer les nouveaux itrs xk+1 et k+1 selon : xk+1 = xk + dx et k+1 = k + d , o est un rel que l'on diminuera jusqu' vrier une condition de descente du genre Wolf ou Armijo (voir [2]). L'existence d'un tel coecient est garantie par l'ingalit (x)dx < 0, qui se traduit par : > 0 tel que

|| < , (x + d) (x) = (x)d + o() < 0

Un choix logique pour la fonction de mrite pendant la Recherche Linaire serait (x), nanmoins, cette fonction ne dpend pas de , alors que le coecient dcide du pas de dplacement pour les deux variables (x, ). Certains auteurs ont donc propos d'ajouter (x) une fonction dpendant de pour la solution de (P ), une possibilit est de choisir comme . Et puisque doit vrier i = ci (x) fonction de mrite 2 : (x, ) = (x) + (x, ) o :
n

(x, ) = T c(x)
k=1

log(k ck (x))

Il est possible de montrer (voir l'article [1]) que (dx , d ) reste une direction de descente pour (x, ) aussi, ce qui justie la Recherche Linaire avec cette fonction de mrite. En utilisant la fonction (x, ) avec le critre de descente d'Armijo, l'tape de la Recherche Linaire peut s'crire : On dispose au dbut d'un pas initial 0 et des constantes , 1 et 2 .
RL :

1. On calcule (x + dx , + d ). 2. Si : (x + dx , + d ) < (x, ) + (x, )T d , on accepte le pas .

2 Le
minimum de la fonction

u u log(u)

est atteint pour

u=

avec

>0

et

> 0.

3. Sinon, on diminue le pas [1 , 2 ] et on revient vers 1 En gnral, on prend : 0 = 1, = 104 , 1 = 0.01 et 2 = 0.99 La plus simple mthode pour diminuer le pas et de choisir 1 = 2 = , ie : j = j . Ceci pose un problme pour le choix de qui, s'il est trop grand, peut imposer un nombre trop important d'itrations pour trouver un bon pas, et s'il est trop petit, peut causer une volution trop lente des itrs et peut-tre mme une stagnation de l'algorithme. Ceci peut tre expliqu en partie par le fait que la diminution de est faite, dans ce cas, sans prendre en compte les informations disponibles sur la valeur des fonctions et de leurs drives depuis le dbut de la RL. En eet, chaque itration j de la RL on dispose de la valeur de (x + i dx , + i d ) pour 0 i j et (x, )T d qui peuvent tre exploites pour la diminution de . Une possibilit donc est de construire un modle quadratique ou cubique, selon le nombre d'informations disponibles, pour la fonction relle () = (x + dx , + d ) qui sera minimiser, par rapport la variable , dans l'intervalle [1 , 2 ] o est la dernire valeur du pas de descente qui ne satisfait pas le critre de descente.

1.2.3 Actualisation de l'information du second ordre

Le calcul de la direction de descente (dx , d ) implique l'valuation du hessien du lagrangien qui n'est pas disponible en gnral sous forme analytique et qui ncessite donc plusieurs valuations de la fonction objectif f pour l'approximer par dirences nies. Il serait donc intressant du point de vue de la vitesse de l'algorithme d'utiliser une mthode qui ne ncessite pas l'valuation des drives secondes. La mthode la plus rpandue est celle de Quasi-Newton qui consiste approximer le hessien de par une matrice M strictement positive qui sera mise jour aprs chaque itration par une formule dite BFGS, le principe tant de n'utiliser que les informations dj calcules ie : (xk+1 , k+1 ), (xk , k ) et Mk . Cette formule est base sur le dveloppement de Taylor avec reste intgral :
1

k =
0

(xk + tk )dt k
x

avec : k = xk+1 xk et k =

(xk+1 , k+1 )

(xk , k+1 ).

Pour approximer le hessien 2 au point xk+1 , une possibilit est donc de choisir une matrice T Mk+1 qui vrie : k = Mk+1 k et qui vrie, comme le hessien, Mk+1 = Mk+1 . Ces deux quations sont insusantes pour dterminer Mk+1 , car on dispose de moins d'quations que d'inconnues. Pour la formule BFGS, on impose en plus qu'elle soit dnie positive et qu'elle minimise une "certaine distance" la matrice, dj calcule, Mk (voir [2]). On arrive alors la formule :

Mk+1 = Mk

T k T Mk k k Mk + T k T k Mk k k k

(F ormule BF GS)

T Il est possible de montrer que Mk+1 est dnie positive si Mk l'est et k k > 0. Il est interessant noter que Mk+1 Mk est un matrice de rang 1, ce qui confre la suite (Mk )k une certaine stabilit aux cours des itrations. Plusieurs changements de la formule BFGS ont t proposes pour rpondre T direntes exigences pratiques. Par exemple, la condition k k > 0 n'est pas garantie dans la majorit des problmes d'optimisation. Heureusement, la correction de Powell donne une solution cet obstacle.

Correction de Powell :

T Pour garder la proprit de la dnie-positivit de Mk+1 mme quand k k 0 (ce qui peut arriver en cas de non convexit du lagrangien), il est souvent fait appel une heuristique appele correction

P de Powell, qui consiste remplacer k par un vecteur : k = k + (1 )Mk k . Le paramtre est T P T pris aussi grand que possible dans [0, 1], de telle sorte que k k k Mk k , pour une constante dans ]0, 1[ (typiquement = 0.2). On trouve : T T si k k k Mk k 1 T Mk k = sinon (1 ) T k T k Mk k k k
Mise l'chelle rpte :

Une deuxime modication de la formule BFGS ci-dessous est d'eectuer une mise l'chelle de la matrice Mk pour mieux estimer le hessien. En eet, la F ormule BF GS est une approximation par une matrice constante du hessien 2 (xk + tk ), or ceci n'est plus judicieux quand les deux points xk et xk+1 sont trs loigns, ce qui peut arriver en cas de mauvais choix du point initial. Il est donc souhaitable de prendre en compte une possible volution rapide des itrs lors l'estimation du hessien. Le choix le plus utilis est :

Mk+1 = k Mk
avec

T M k k k Mk TM k k k T k k TM k k k

T k k T k k

k =

Ce coecient a t modi par Al-Baali (voir l'article [11]) : k = min

rsultats numriques obtenus dans [11] avec ce choix montre qu'il est plus performant.

T k k , TM k k k

et les

1.2.4 Algorithmes A et A

Nous pouvons enn prsenter l'algorithme en pseudo-code. Comme c'est expliqu prcdemment, le problme initial P (algorithme global A) est une suite de problmes P (algorithme A ) o l'on diminuera le paramtre au fur et mesure.

Remarque :
La mthode la plus simple pour diminuer le paramtre , est : + = , avec > 1 pour que la suite k converge vers 0. Des techniques de mise jour plus sophistiques existent, mais pour notre implmentation nous avons choisi un simple critre, utilis dans d'autres algorithmes reconnus : Si le problme (P ) a t rsolu en moins de trois itrations, on prend : + = 100 , sinon + = 5

Algorithm 1 while

Algorithme (A ) || f (x) c(x)|| > Or ||C(x) e|| >

do

1) Calcul de la direction de descente :

d = (dx , d ) solution de : (M + c(x)C(x)1 c(x)T )dx = f (x) + c(x)C(x)1 e c(x)T dx + C(x)d = e C(x)
2) Recherche Linaire avec critre d'Armijo : while

(x + dx , + d ) > (x, ) + (x, )T d [1 , 2 ]

do

end while

x+ = x + dx + = + d
3) Actualisation de l'information du second ordre : Formule BFGS

= x+ x = x (x+ , + ) x (x, + ), avec correction de Powell si necessaire. T k k , 1 k = min T k Mk k T Mk k k Mk k T Mk+1 = k Mk + T k TM k k k k k

end while

Algorithm 2

Algorithme global A On dispose d'une valeur initial 0 (On peut simplement prendre 0 = 1) et d'un point initial (x, ) strictement admissible.
while

|| f (x )

c(x ) || > tolerance Or ||C(x ) || > tolerance

do

Diminuer le paramtre .

(x , ) sera le point initial du nouveau problme P+ (x+ , + ) solution de P+ (sortie de l'algorithme A+ )

end while

1.2.5 Convergence de l'algorithme global A

Le thorme suivant assure la convergence da l'algorithme bonne solution.
Thorme 1.2.1 A,

sous certaines conditions vers la

Convergence de l'algorithme global A (voir [1]). Supposons que :  f et c1 , ..., cm sont convexes et qu'au moins une de ces fonctions est fortement convexe.  f et c sont C 1,1 .  x Rn tel que c(x) > 0. Alors l'algorithme A gnre une suite borne {zi } et sa limite est une solution du problme primal dual (P).
Une preuve peut tre trouve dans l'article [1].

1.3 Rsultats Numriques

Nous prsenterons dans cette partie plusieurs exemples de problmes d'optimisation sur lesquels on testera l'algorithme A. Ils sont des cas classiques auxquels on soumet tout algorithme d'optimisation, car ils prsentent des particularits qui rendent leur minimisation dicile. Ils sont soit issus d'autres branches comme la mcanique ou la logistique, soit des problmes articiels dans le sens o ils ont t conus pour tester la robustesse des algorithmes. Nous allons leur appliquer l'algorithme A en utilisant plusieurs points initiaux, puis estimer la meilleure solution du problme trouve, la moyenne des solutions et la pire. Ceci nous permettra d'tudier la convergence de l'algorithme et sa dpendance au point initial choisi, le but nal tant de pouvoir utiliser l'algorithme pour la calibration du modle d'Heston. Pour commencer, nous avons test l'algorithme A avec un ensemble de problmes tests G1, G4, G6 , G7, G8, G9, G10 (voir l'article [12]), qui sont rappels dans l'annexe B. Les caractristiques de ces problmes tests sont susamment diverses pour couvrir plusieurs sortes de dicults que rencontrent les algorithmes d'optimisation. D'autres problmes sont tudis dans la section suivante. Nous avons test notre algorithme en utilisant 20 points initiaux dirents et rapport, dans la table en bas, la meilleure solution connue dans la littrature3 , la meilleure solution trouve par l'algorithme, la moyenne et la pire solution, en plus de l'cart type des rsultats :

Test Problem G1 G4 G6 G7 G8 G9 G10

Best Known -15 -30665.539 -6961.81388 24.3062091 -0.095825 680.6300573 7049.3307

Best -14,999997 -30665,538672 -6961.813876 24.30620906 -0,095825 680.6300573 7049.248020769

Worst -14,851949 -30665,535061 -6961.813875 24.30621290 -0,095825 680.6300893 7049.248500528

Average -14,980565 -30665,538343 -6961.813875 24.30620928 -0,095825 680.6300610 7049.248091997

Std Deviation 0,039242 0,001089 0.00001 0.00000098 0 0.0000000 0.000214125

Les rsultats de l'algorithme A sont trs satisfaisant ( part pour le problme G6), et les valeurs minimales sont toutes retrouves pour tous les problmes. L'algorithme arrive mme amliorer la valeur minimale pour le problme G10. Pour ce problme, la meilleure valeur minimale connue est 7049.3307, mais nous trouvons 7049.2480, et ce pour :

x = (579.30668, 1359.97067, 5109.97067, 182.017610, 295.60117, 217.98230, 286.41653, 395.60117)

Dans la suite nous tudierons d'autres problmes d'optimisation issus des dirents domaines de l'ingnierie.
3 Valeurs
reprises de l'article [8], table 2.

10

1.3.1 Fonction de Rosenbrock

Cette fonction est un test classique auquel on soumet les algorithmes d'optimisation sans contraintes.
N 1

Elle est dnie sur RN par : f (x) =

i=1

100(xi+1 x2 )2 + (1 xi )2 . Elle admet un minimum global i

unique au point x avec x = 1 et f (x ) = 0. En dimension 2, le graphe de la fonction est caractris i par une valle autour de x , ce qui peut induire un algorithme en stagnation.

Fig.

1.1  Graphe de la fonction de Rosenbrock

Nous avons test notre algorithme avec cette fonction, sous les contraintes de bornes xi [2, 2], pour 100 points initiaux choisi alatoirement, et pour direntes valeur de N . Dimension N 2 5 10 20 30 Pourcentage d'checs 0% 9% 15% 17% 20%

1.3.2 Welded Beam Design Problem

La fonction objectif pour ce problme a quatre variables x = (x1 , x2 , x3 , x4 ) et peut tre exprim sous la forme (voir l'article [9]) :

11

 Minimiser f (x0 , x1 , x2 , x3 ) = 1.10471x2 x2 + 0.04811x3 x4 (14.0 + x2 ) 1 Sous les contraintes : g1 = t tmax 0 g2 = s smax 0 g3 = x1 x4 0 g4 = 0.10471x2 + 0.04811x3 x4 (14.0 + x2 ) 5.0 0 1 g5 = d dm ax 0 g6 = P Pc 0 avec : P = 6000; L = 14; E = 30e6; G = 12e6 tmax = 13600; smax = 30000; dmax = 0.25 M = P (L + x2 ); R = 0.25(x2 + (x1 + x3 )2 ) 2 22 x2 2 J = 2 x1 x2 ( 12 + 0.25(x1 + x3 )2 ) E Pc = 4.013E x3 x3 (1 0.25x3 LG ) 2 4 6L P t1 = 2x x ; t2 = M R/J 1 2 t = t2 + t1 t2 x2 /R + t2 2 1 6P L s = x4 x2 3 4P L3 d = Ex4 x3
3

(1.2)

En lanant l'algorithme A avec 30 points initiaux, nous obtenons les statistiques suivantes : BEST SOLUTION WORST SOLUTION MEAN STANDARD DEVIATION 1.724852309 1.724862549 1.724853257 2.44719E-06

Nous trouvons donc la valeur minimale suivante pour f : 1.724852309, alors que la meilleure valeur trouve dans [9] est 1.728226. Nous trouvons donc une meilleure valeur et cela pour les variables :

x = (0.205729627, 3.470489088, 9.03662404, 0.205729645)

1.3.3 Pressure Vessel Design Problem

La fonction objectif pour ce problme est le cot de fabrication d'un cylindre contenant de l'air sous pression de 3000 psi (voir gure 1.2), les contraintes reprsentant des limitations sur ses dimensions.

2 2 2 Minimiser f (x0 , x1 , x2 , x3 ) = 0.6224x0 x2 x3 + 1.7781x1 x2 + 3.1661x0 x3 + 19.84x0 x2 Sous les contraintes : x1 0.00954x2 0 x0 0.0193x2 0 (x2 x3 + 4 x3 ) 1296000 0 2 2 3 240 x3 0

(1.3)

La valeur minimale de la fonction objectif trouve dans l'article [8] est de fmin = 5868.764836. Comme dans beaucoup de problmes d'optimisation, la convergence de l'algorithme dpend du point initial cause de l'existence ventuelle de minima locaux. Nous avons essay notre algorithme avec dirents points initiaux, en variant x2 et x3 et en xant x0 et x1 . 12

Fig.

1.2 

Nous obtenons pour x0 = 1.5 et x1 = 0.9 les rsultats donns dans la gure 1.3.

Fig.

1.3  Rsultats pour le Vessel Design Problem

Et en variant aussi les variables x0 et x1 , nous obtenons pour 100 points initiaux les statistiques suivantes : BEST SOLUTION WORST SOLUTION MEAN STANDARD DEVIATION 5804,3762 6640,1494 5862,5904 165,3692

Nous obtenons donc une valeur minimale plus petite que celle trouve dans l'article [8], et cela pour x = (0.727590929, 0.359648573, 37.69901189, 239.99999997). L'ensemble des tests auxquels on a soumis l'algorithme A donne des rsultats assez prometteurs pour son utilisation en calibration du modle d'Heston. Comme cette calibration se fait souvent sur 13

les prix des options vanilles, il est essentiel de disposer d'une mthode d'valuation de ces instruments qui soit stable numriquement. Dans ce sens, nous allons, dans la prochaine partie, rappeler le cadre du modle d'Heston et la formule semi explicite pour l'valuation d'une option europenne. Ensuite, nous discuterons des problmes numriques rencontrs lors de son utilisation et les solutions apportes ces problmes.

14

Chapitre 2
Modle d'Heston : Evaluation

2.1 Prsentation du modle

Le modle d'Heston a t introduit pour rpondre plusieurs constatations exprimentables :  L'volution alatoire de la volatilit a une partie dcorrle de la dynamique du sous-jacent, ce qui justie la modlisation de la volatilit par un processus stochastique part entire, et pas seulement comme fonction dterministe du sous-jacent, comme c'est le cas dans les modles volatilit locale.  Les phnomnes de smile remarqus dans le march et la dynamique temporelle de la surface de volatilit implicite, que le modle de volatilit locale ne permet pas de retrouver, contrairement au modle d'Heston. La dynamique propose par Heston [3] a l'avantage d'expliciter cette corrlation avec un seul paramtre et de donner une volution raliste la volatilit. Il modlise la volatilit comme un processus de diusion vriant une EDS de type Cox-Ingersoll-Ross. La dynamique du modle est de la forme : dSt = St rdt + St Vt dWt dVt = ( Vt )dt + Vt dZt d < Wt , Zt >= dt avec valeurs initiales S et V0 . Elle est determine par 5 paramtres (, , , V0 , ).

(2.1)

2.1.1 Proprits thoriques

L'existence du processus Vt ne peut tre assure par les thormes classiques (bass sur la mthode de Picard), car la racine carre n'est pas lipchitzienne en 0. Par contre on peut utiliser un thorme de Yamada-Watanabe qui nous permet d'obtenir l'existence et l'unicit de ce processus (Voir annexe A). Vs t t Ensuite, le processus St est donn par : St = S0 e 0 (r 2 )ds+ 0 Vs dWs
Moments du processus de la variance :

En plus, on peut montrer 1 que, grce la croissance sous-linaire des coecients des l'EDS, on a un contrle uniforme sur les moments de Vt et de Xt = log(St ), pour tout p > 0 :

E
1 voir

Sup |Xt |p + |Vt |p C(p, T )(1 + |X0 |p + |V0 |p )

0tT

l'annexe A pour une preuve

15

Ceci nous permet par exemple de calculer les moments de Vt .

t t En eet : Vt = V0 + 0 ( Vs )ds + 0 Vs dZs , et d'aprs le contrle sur les moments de Vt , t t l'intgrale 0 Vs dZs est une martingale, donc de moyenne nulle, d'o E[Vt ] = E[V0 ]+ 0 (E[Vs ])ds. La moyenne vrie donc une quation direntielle simple rsoudre dont la solution est :

E[Vt ] = + (V0 )et

De manire similaire on trouve que :

V ar(Vt ) E[ln St ]

2 V0 2 t e (1 et ) + (1 et ) 2 1 1 et = ln F t (V0 ) 2 2

avec F = erT S est la valeur forward de S .  La moyenne de Vt tends vers quand t tends vers l'inni, avec une vitesse controle par , ce qui justie le nom de moyenne de la variance pour et taux de retour la moyenne pour .  La variance de Vt augmente avec la volatilit de la variance et diminue avec . On signale enn qu'il est possible de montrer que le processus V reste strictement positif, si la condition 2 2 > 0 est vrie. Cette ingalit est appele condition de Feller.
Fonction caractristique du log-spot :

La loi de ST est connue travers la transforme de Fourier de ln ST : (u) = E[eiu

ln ST

].

En cherchant (u) sous la forme (u) = eC(T,u)+D(T,u)V0 +iu ln F , et en utilisant l'EDP vrie par (u) par rapport aux variables S0 , V0 et T , on montre que C(T, u) et D(T, u) vrient deux quations de Riccati faciles rsoudre. On trouve nalement :

C(T, u) =

( ui + d(u))T 2 ln ui + d(u) 2 et d(T, u) =

c(u)ed(u)T 1 c(u) 1

D(T, u) =
avec :

ed(u)T 1 c(u)ed(u)T 1 (ui )2 + iu 2 + 2 u2

c(T, u) =

ui + d(u) ui d(u)

Vu que la racine carre d'un nombre complexe est une fonction multiforme, il nous faut, pour calculer z , choisir une forme qui sera forcement discontinue sur une demi droite. On choisit la forme principale (qui admet comme demi-droite de discontinuit R ), et en plus on peut choisir soit z soit z . Le choix sera motiv par la ncessit du calcul du logarithme complexe, qui est source de dicult numrique comme c'est expliqu dans la section suivante.

2.2 Prix d'un Call : Formule semi-explicite

Connaissant la transforme de Fourier du log-spot, le calcul du prix d'un call peut se faire facilement en utilisant la formule de Fourier inverse pour retrouver la densit du log-spot.

16

En eet, si l'on nomme p(.) cette densit, alors : p(s) =

1 2

+ ius e (u)du.

D'autre part : C(T, S, K) = E[erT (ST K)+ ] = erT E[ST 1ST K ] + erT K E[1ST K ] or :
+

E[1ST K ]

p(s)1sln K ds 1 2 1 2
+ + +

= =

eiu s (u)1sln K dsdu

+

eiu

ln K

(u)du

eiux 1x0 dx

et en utilisant la transforme de Fourier inverse de la fonction de Heaviside :

+

eiux 1x0 dx = (u) +

1 iu f (x)(x)dx = f (0) on trouve :

, o est la distribution de Dirac, et avec la proprit

E[1ST K ]

= =

1 1 + 2 2 1 1 + 2
0

+ +

eiu ln K (u) du iu eiu ln K (u) iu du

dP ST = , on obtient : E[ST 1ST K ] = E[1ST K ], et on dP F peut alors utiliser la formule trouve ci-dessus, sachant que la fonction caractristique de ST sous P 1 est gale : (u) = F (u i). On trouve alors :
En utilisant le changement de probabilit

E[ST 1ST K ] =
Donc :

1 F + 2

+ 0

eiu ln K (u i) iu

du

C(T, S, K) = erT
avec :

F K 1 + 2 eiu

(F f1 (u) K f2 (u))du
0 ln K

f1 (u) =

(u i) iuF
ln K

f2 (u) =

eiu

(u)

iu

En pratique, le calcul de C(T, S, K) soulve deux dicults : La premire est le calcul de la fonction caractristique, qui fait intervenir le logarithme complexe, chose qui n'est pas vidente et qui peut induire des instabilits numriques. Le second point est le choix d'un schma d'intgration numrique qui soit robuste et adapt aux fonctions en jeu.

2.2.1 Instabilit numriques : Logarithme complexe

Le problme de l'valuation du logarithme complexe se fait remarquer exprimentalement en traant la fonction f (u) = F f1 (u) K f2 (u) (voir la gure 2.1). Thoriquement, cette fonction est continue sur ]0, [ (sa continuit en 0 est dmontre en annexe A), par contre son graphe rvle des 17

discontinuits pour certains paramtres. En tudiant de prs la formule d'valuation d'un call, on remarque que l'valuation du prix d'un call passe par l'intgration de f (u) sur [0, +[, ce qui ncessite le calcul de log
1cedT 1c

Le problme vient du fait que le logarithme complexe est une fonction multiforme : elle admet plusieurs dterminations. En eet, si w C\{0}, alors z = log(w) = x + iy est dtermin par ez = w. ie : ex = |w|, soit x = log |w| et d'autre part :

y = arg(w)

Il existe plusieurs dterminations, mais elles dirent les unes des autres d'une multiple entier de 2i . On montre que la fonction arg n'admet pas de dtermination continue sur C\{0}. On est oblig d'en choisir une, ce qui correspond choisir une demi droite travers laquelle arg est discontinue. En gnral le choix est x sur la demi droite R et la discontinuit arrive quand Im(z) change de signe avec (z) < 0. 1 cedT traverse la demi droite R quand u parcourt Il est donc ncessaire d'tudier si G(u) = 1c [0, [. Ceci dpend du choix du signe de d(T, u) = (ui )2 + iu 2 + 2 u2 , o la racine est calcule avec sa forme principale (qui admet R comme axe de discontinuit). En eet, l'article [5] montre qu'avec le choix d(T, u) = d1 (T, u) = + (ui )2 + iu 2 + 2 u2 , G traverse l'axe R pour des maturits assez grandes, c'est dire que log(G) calcul avec sa dtermination principale subit une discontinuit. Une solution ce problme, propose par l'article [4], consiste corriger la discontinuit due la partie arg du logarithme en ajoutant des multiples de 2 pour rendre la fonction continue. Par contre, le seul choix d(T, u) = d2 (T, u) = (ui )2 + iu 2 + 2 u2 assure la continuit de log(G(u)). En eet, le thorme 3 de l'article [5] montre qu'avec le choix prcdent, G(u) et G(ui) ne traverse pas l'axe des rels ngatifs pour tout u R+ pour n'importe quel jeu de paramtres. Ces problmes numriques se manifestent sous forme de discontinuit de la fonction f (u) que l'on doit intgrer sur [0, [, comme le montre la gure en bas o l'on a trac f (u), en utilisant la forme d1 (T, u), sans la correction propose dans l'article [4], puis avec le choix d2 (T, u) qui ne ncessite pas de correction.

Fig.

2.1  Discontinuit de au logarithme complexe.

18

2.2.2 Intgration Numrique

Le deuxime point sensible dans l'estimation de C(T, S, K) est l'utilisation d'un schma de discr+

tisation numrique pour calculer l'intgrale

b 0

(F f1 (u) K f2 (u))du.

En gnral, pour estimer

m a

f (x) dx, on a recours des formules d'approximation du type :

I=
i=0

ai f (xi ) dans lesquelles les noeuds xi appartiennent [a, b] et les poids ai sont des rels.

Un exemple de schma simple est la mthode des trapzes. Elle consiste interpoler la fonction sur chaque subdivision par un polynme de degr 1 :
b a n (f ),

ba f (x) dx = n

f (a) + f (b) + 2
3

n1

f
k=1

a+k

ba n

n (f )

est l'erreur de quadrature et vaut : (ba) f () pour un [a, b] 12n2

Mais comme c'est signal dans l'article [4], la forme des fonctions f1 et f2 peut trs bien tre de nature exponentielle, ou bien trs oscillatoire, selon le choix des paramtres. Il est donc essentiel de choisir un schma d'intgration plus avanc. Pour cela, nous utiliserons le schma adaptive Gauss-Lobatto propos dans l'article [6]. Il bas sur l'valuation itrative de l'intgrale sur un grille adaptative. En eet, l'intervalle [a, b] est rcursivement divis en sous intervalles sur lesquels, l'intgrale est value par une des rgles de quadrature.

19

Chapitre 3
Modle d'Heston : Calibration

La pertinence d'un modle est souvent value par plusieurs critres parmi lesquels la capacit des prix des options vanilles gnrs par le modle coller aux prix du march. En eet, le but de la calibration est de pouvoir utiliser le modle, une fois calibr, pour calculer le prix des options exotiques et la stratgie de leur couverture. Or, les options vanilles sont souvent utilises comme des outils de couverture dans ces stratgies, d'o l'intrt d'avoir un modle qui colle assez bien ces options-l. C'est d'ailleurs l'une des raisons pour laquelle les modles volatilit stochastique ont t introduits aprs les modles volatilit locale. En eet, comme nous l'avons soulign au dbut du chapitre prcdent, la dynamique temporelle du smile prdite par les modles volatilit locale ne correspond pas celle observe dans le march. Le principe donc est, faute de pouvoir coller parfaitement l'ensemble des prix rels, de retrouver des paramtres qui minimise la distance entre les prix du modle et ceux du march, ce qui revient un problme d'optimisation que l'on peut crire, avec la notation x = (, , , V0 , ) :

min f (x) xX
o f est une fonction d'erreur et X un ensemble dans lequel on cherche nos paramtres optimaux. Comme dans tout problme d'optimisation, la premire tape est d'tudier la dpendance de la fonction f par rapport la variable x. On va donc dans la partie suivante tudier l'inuence des dirents paramtres sur la nature des courbes des prix.

3.1 Volatilit Implicite : Eet des paramtres

Avant de discuter de la calibration du modle, il est intressant d'tudier le rle de chacun des paramtres dans la construction des nappes de volatilits implicites, outil essentiel utilis par les traders pour analyser la structure en terme des prix. Une ventuelle corrlation entre les paramtres, par exemple, peut permettre de rduire leur nombre pendant la calibration, en en xant certains.  Le premier paramtre auquel on s'intresse est la vitesse de retour la moyenne . D'aprs le rle de dans la dynamique de la volatilit, il est logique d'tudier sparment son eet pour des courtes et longues maturits. En eet, l'eet de est plus important pour les grandes maturits. Comme le montre la gure 3.1, le paramtre joue sur la convexit du smile : l'augmentation de diminue la convexit de la nappe de volatilits. Nous remarquons aussi que ce changement est plus prononc pour les volatilits ATM.  La volatilit de la variance a un eet similaire celui de : la gure 3.2 montre que contrle lui aussi la convexit du smile.

20

Fig.

3.1  Eet de la vitesse de retour la moyenne sur la convexit de la volatilit implicite

Fig.

3.2  Eet de la volatilit de la variance sur la convexit de la volatilit implicite

 Le troisime paramtre est la volatilit long terme : intuitivement, une augmentation de remonte le prix, similairement au paramtre dans le modle de Black-Scholes. Ceci se rete dans la gure 3.3, qui montre que contrle le niveau du smile. Le mme rle est jou par la variance initiale V0 .  Le dernier paramtre est la corrlation : son changement entrane un changement de la symtrie 21

du smile, et a un eet oppos sur les prix OTM et ITM. En eet, augmenter la corrlation augmente les prix des Calls OTM et diminue les prix des Calls ITM. Ceci rete l'eet de la corrlation sur la skewness de la distribution des log-rendements, qui a le mme signe que la corrlation. Ces direntes remarques peuvent nous suggrer de considrer et comme une mme variable, de mme pour et V0 , ce qui diminue le nombre de variables calibrer 3 paramtres, ce qui peut tre intressant pour rduire le temps de calcul lors de l'optimisation. Mais pour plus exibilit, nous ne ferons pas cette hypothse et nous conserverons les 5 paramtres spars pendant la calibration

Fig.

3.3  Eet de la volatilit long terme sur le niveau de la volatilit implicite

Fig.

3.4  Eet de la corrlation sur le skew de la volatilit implicite

22

3.2 Calibration : Introduction

Comme nous l'avons soulign dans la partie prcdente, le principe de la calibration est de minimiser la distance entre les prix du modle et ceux du march, ce qui revient un problme d'optimisation que l'on peut crire, avec la notation x = (, , , V0 , ) :

min f (x) xX

(3.1)

o f est une fonction d'erreur et X un ensemble qui reprsente les contraintes sur les paramtres. En l'occurrence, les paramtres doivent vrier les contraintes de bornes suivantes : > 0, > 0, > 0, V0 > 0 et 1 < < 1. En plus, nous imposerons la condition de Feller 2 2 > 0 comme contrainte, ce l'instar de plusieurs recherches qui ont montr que cette condition est ncessaire pour la stabilit des prix des options exotiques. Il existe plusieurs choix possibles pour la fonction d'erreur f , qui rpondent dirents impratifs. Le choix le plus utilis est la fonction d'erreur quadratique que l'on notera rmse, pour root mean square error 1 :

rmse(x) =

(M arket price M odel price)2 number of options options

Clairement, cette fonction donne plus de poids aux options chres, c'est--dire les options ITM et grande maturit, au dtriment des options trs OTM et courte maturit. Une solution serait de mettre des poids aux prix, ou bien d'utiliser d'autres fonctions d'erreurs qui value la dirence relative entre le prix du march et celui du modle, comme l'average relative percentage error :

arpe(x) =

1 |M arket price M odel price| number of options options M arket price

Il est aussi possible d'utiliser la fonction average absolute error dnie par :

aae(x) =

1 |M arket price M odel price| number of options options

Le problme que nous devons rsoudre s'crit alors, en utilisant par exemple la fonction rmse : min rmse(, , , V0 , ) > 0, > 0, > 0, V0 > 0 et 1 < < 1 (3.2) 2 2 > 0 Dans la suite, nous utiliserons les trois fonctions pour direntes sortes de donnes. En premier lieu, nous tudions la capacit du modle retrouver les paramtres avec lesquels une nappe de prix a t gnre. Enn, nous utiliserons l'algorithme march.
A

pour calibrer le modle d'Heston sur des donnes relles du

3.2.1 Application aux prix gnrs par le modle d'Heston

Avant d'utiliser l'algorithme A pour la calibration du modle d'Heston sur des donnes relles du march, on soumet l'algorithme un dernier test qui consiste retrouver les paramtres avec lesquels une nappe de prix a t gnre.
1 On
suit dans cela les notations de l'article [7]

23

Avec les paramtres (, , , V0 , ) = (1.5, 0.5, 0.16, 0.16, 0.1), nous avons calcul un ensemble de prix avec 6 maturits de 0.5 Y 15 Y , chacune avec 7 strikes de 40% 150%, en xant S0 = 100, r = 0.03. Pour tudier la dpendance de l'algorithme A au point initial, nous avons utilis 10 dirents jeux de paramtres initiaux. Nous donnons dans le tableau en bas, l'cart type des paramtres trouvs aprs la calibration ralise avec la fonction rmse :

rmse Ecart-type
Valeurs d'origine Valeurs trouves

5.77E-03 1.5 1.8238

1.18E-03 0.8 0.7575

1.50E-05

V0 1.97E-05 -0.1 -0.1432

2,25E-04
Cond.Feller -0.16 1.69E-4

0.16 0.1573

V0 0.16 0.1611

On voit clairement que nous arrivons retrouver les paramtres originaux, avec une bonne stabilit par rapport au point initial choisi. On signale aussi que l'on obtient une meilleure stabilit avec la fonction arpe.

arpe Ecart-type

4.05E-03

4.45E-04

1.84E-05

V0 9.54E-06

7.74E-05

Par exemple, avec le point initial suivant : (, , , V0 , ) = (4, 1, 0.4, 0.4, 0.5), nous obtenons les paramtres suivants : (, , , V0 , ) = (1.49967014, 0.50012982, 0.16000544, 0.16001793, 0.10013121).

3.2.2 Calibration du modle d'Heston sur des donnes de march

D'aprs la partie prcdente, nous avons montr que l'algorithme donne des rsultats prometteurs pour son utilisation pour la calibration du modle d'Heston sur des prix du march. Nous allons donc considrer dans cette partie des prix rels et, pour des besoins de comparaison, nous utiliserons les donnes de l'article [7], composes de 144 Call europens sur l'Eurostoxx 50 index le 07/10/2003, dont la valeur est 2461.44, et pour des maturits de 0, 0361 Y 5, 1639 Y . Nous prenons la mme valeur que dans l'article pour le taux d'intrt sans risque r = 0.03. (voir la table 4 de l'article [7]) Comme pour la partie prcdente, nous tudierons la qualit de la calibration et sa dpendance du choix du point initial. En eet, nous avons lanc notre algorithme avec 10 dirents paramtres initiaux en utilisant les trois fonctions d'erreur rmse, arpe et aae. Dans le tableau en bas, nous donnons la valeur minimale de la fonction d'erreur obtenue aprs la calibration et, pour comparer, nous donnons aussi la valeur trouve dans l'article [7] : Article [7] Nos rsultats

rmse 3.0281 2.8365

arpe 0.0174 0.0101

aae 2.4264 2,3209

Nous remarquons que l'algorithme A trouve une meilleure valeur pour les trois fonctions d'erreur par rapport l'article [7]. La gure 3.5 montre que le modle d'Heston arrive reconstituer la nappe des prix avec une bonne qualit. Le signe plus reprsente les prix donns par le modle et les cercles ceux du march. En plus, l'algorithme est stable par rapport au point initial comme le montre le tableau suivant qui donne l'cart type des paramtres trouvs avec les 10 points initiaux : Ecart-type

1.68E-04

7.64E-06

4.91E-06

V0 1.49E-06

6.05E-05

Nous donnons aussi les paramtres donnant la meilleure valeur pour les trois fonctions d'erreur :

24

Fig.

3.5  Calibration du modle d'Heston

rmse arpe aae

0.49732481 0.473203521 0.51500004

0.265879774 0.255870681 0.265000991

0.071072575 0.069805804 0.070417831

V0 0.06438063 0.063142791 0.06454543

-0.733659533 -0.692745768 -0.748583081

Comme nous l'avons soulign dans l'introduction, l'utilisation de la fonction d'erreur quadratique donne plus de poids aux options chres, ie de longue maturit ou trs dans la monnaie, par rapport aux options courte maturit et trs hors de la monnaie. Ceci implique une mauvaise calibration sur ces |M odel price M arket price| dernires options, comme le montre la gure 3.6 qui trace l'erreur relative M arket price aprs la calibration. L'erreur relative est plus grande pour les grands strike (options trs OTM).

Fig.

3.6  Erreur relative pour rmse

Ce problme peut tre rsolu en partie en utilisant la fonction d'erreur relative arpe dans la calibration. La gure 3.7 nous montre que l'erreur relative sur les options OTM a t presque divise par deux.

25

Fig.

3.7  Erreur relative pour arpe

Conclusion
Le sujet de ce rapport a t l'utilisation d'un algorithme d'optimisation pour la calibration du modle d'Heston. Nous avons donc tudi un algorithme d'optimisation qui utilise la mthode des Points Intrieurs et les techniques de quasi-Newton. Les dirents tests auxquels on a soumis l'algorithme ont montr sa robustesse par rapport au choix du point de dpart. L'utilisation de cet algorithme pour la calibration du modle d'Heston a donn des rsultats satisfaisants, que ce soit pour qualit de la convergence ou pour la stabilit des rsultats de la calibration par rapport aux paramtres initiaux. Par contre, la qualit de la calibration est relativement limite pour les courtes maturits, ce qui est un des inconvnients du modle d'Heston. Une possibilit donc est d'utiliser un modle plus avanc, comme celui de Bates ou avec un processus de Lvy gnral.

26

Bibliographie

[1] P. Armand, J. C. Gilbert ET S. Jan-Jgou (2000). A feasible BFGS interior point algorithm for solving strongly convex minimization problems, SIAM J. Optim. 11, no. 1, 199-222. [2] J.Ch. Gilbert.(2007) Optimisation Direntiable - Thorie et Algorithmes, Syllabus de cours l'ENSTA. http ://www-rocq.inria.fr/ gilbert/ensta/optim.html. [3] Heston, S.L.(1993) A closed-form solution for options with stochastic volatility with applications to bond and currency options, Review of Financial Studies , vol. 6, no. 2, pp. 327-343. [4] Kahl, C. and Jackel, P. (2005) Not-so-complex logarithms in the Heston model, WilmottMagazine, September 2005. [5] Albrecher, H., Mayer, P., Schoutens, W., Tistaert, J. (2007) The little Heston trap, Wilmott Magazine, January Issue (2007) 83-92. [6] Gander, G. and Gautschi, W. 2000 Adaptive quadrature-revisited, BIT Numerical Mathematics 40, 1, 84-101. [7] Schoutens, W., Simons, E. and Tistaert, J. (2004) A Perfect calibration ! Now what ?, Wilmott Magazine, March 2004. [8] Abdel-Rahman Hedar and Masao Fukushima Derivative-Free Filter Simulated Annealing Method for Constrained Continuous Global Optimization [9] C. A. Coello Coello and E. M. Montes(2002) Constraint-handling in genetic algorithms through the use of dominance-based tournament selection Advanced Engineering Informatics, 16 , 193-203. [10] Chandra Sekhar Pedamallu, Linet zdamar, Tibor Csendes An Interval Partitioning Approach for Continuous Constrained Optimization [11] M. Al-Baali (1998) Global and superlinear convergence of a restricted class of self-scaling methods with inexact line searches, for convex functions, Comput. Optim. Appl. 9, no. 2, 191-203 [12] W. Hock and K. Schittkowski Test Examples for Nonlinear Programming Codes, Springer-Verlag Berlin Heidelberg, 1981

27

Annexe A
Modle d'Heston : Rsultats thoriques

A.1 Existence du processus (Vt)t0

L'existence du processus de la variance est assure par le thorme de Yamada-Watanabe qui amliore notablement le thorme classique sur l'EDS. En eet, il n'exige pas que le coecient de la diusion soit localement lipchitzien, comme c'est le cas pour l'EDS du processus CIR :

dVt = ( Vt )dt +

Vt dZt

La fonction . n'tant pas lipchitzienne en 0, le thorme classique ne peut pas s'appliquer, mais on peut utiliser le thorme suivant :
Thorme A.1.1 Thorme de Yamada-Watanabe On considre l'EDS une dimension suivante :

(E)

dXt = b(t, Xt )dt + (t, Xt )dBt X0 = x R

On suppose que :  b et sont croissance linaire, ie :

|b(t, x) + (t, x)| C(1 + |x|), t [0, T ], x R

 b vrie "la condition de Lipchitz locale" ie :

n 1 Kn > 0 tel que t [0, T ], |x| < n, |y| < n |b(t, x) b(t, y)| Kn |x y|

 vrie :
|(t, x) (t, y)| (|x y|), t [0, T ], x R, y R

o :]0, +[]0, +[ est une fonction mesurable telle que :

dz
0

2 (z)

= +

pour tout > 0 Alors, l'quation (E) admet une solution forte unique. L'application de ce thorme avec le choix (x) = x (sachant que | x y| que l'quation vrie par la variance admet une unique solution forte.
28

|x y|) montre

A.2 Contrle uniforme sur les moments de X et V

Comme il a t nonc dans le chapitre 2, les moments de X = log(S) et V sont contrls uniformment dans le sens : Pour tout p 1 :

Sup |Xt |p + |Vt |p C(p, T )(1 + |X0 |p + |V0 |p )

0tT

Ce rsultat est vrai pour tout processus Y dni par une EDS coecients croissance linaire, comme X et V . On considre l'EDS : dYt = b(t, Yt )dt + (t, Yt )dBt (E) Y0 = y R
Thorme A.2.1

(A.1.1). Alors : si T > 0,

On suppose que b et vrient les conditions du thorme de Yamada-Watanabe

p 1 C(p, T ) > 0 tel que : E

Sup |Yt |p C(p, T )(1 + |y|p )

0tT

Preuve : Premier cas : p 2. Pour travailler avec des processus borns, on introduit le temps d'arrt suivant : pour n N : n = inf{t [0, T ], |Yt | > n}.
On a :
t t p

Ytp =
p

y +
0 p1 p

b(s, Ys )ds +
0 p p

(s, Ys )dBs

, et avec l'ingalit : (a + b + c) 3

(a + b + c ), on obtient :
tn p tn p

t T, |Ytn |p 3p1
Donc :

|y|p +
0

b(s, Ys )ds

+
0

(s, Ys )dBs

T n

tn

sup |Yt |
tT n

3
tn 0

p1

|y| + E
0

|b(s, Ys )| ds

+ E sup
tT 0

(s, Ys )dBs

Sachant que Gundy :

(s, Ys )dBs est une martingale, on peut utiliser l'ingalit de Burkholder-Davisp

tn

E
0

T n

p/2

E sup
tT 0

(s, Ys )dBs

(s, Ys )2 ds

On utilise aussi l'ingalit de Holder (la condition p 2 le permet) pour obtenir :

E
0

T n

p/2

CE
0

T n

(s, Ys )2 ds
T n

(s, Ys )p ds
T n

E
0

|b(s, Ys )| ds

CE
0

|b(s, Ys )| ds

29

Par suite, puisque b et sont croissance linaire :

T n

E
0

|b(s, Ys )| ds K

T n

1+E
0

|Ys |p ds

1+E
0

sup |Ys |p ds
sT n

La mme ingalit est valable pour le coecient . On obtient alors :

T

sup |Yt |p A 1 + |y|p +

tT n 0

sup |Yt |p du
tun

o la constante A ne dpend pas de n. Le lemme de Gronwall donne alors :

n0 E

sup |Yt |p A(1 + |y|p )

tT n

Finalement, le lemme de Fatou permet d'avoir :

E sup |Yt |p A(1 + |y|p )

tT

Deuxime cas : si 1 p < 2 On se ramne au cas prcdent en appliquant Holder (remarquer que 2p 2) :
1/2

E sup |Yt |p E sup |Yt |2p

tT tT

1 + |x|2p

C(1 + |x|p )

A.3 Continuit de f1 et f2 en 0
L'tude des fonctions f1 et f2 est indispensable pour mieux cerner les problmes du calcul nu+

mrique de
0

(F f1 (u) K f2 (u))du. Ici nous allons montrer leur continuit en 0. Ceci peut

tre prouv en utilisant l'expression explicite de C et D dans la fonction caractristique (u) = eC(T,u)+D(T,u)V0 +iu ln F (voir [4]). Mais il est plus facile d'utiliser le fait que (u) est la fonction caractristique d'un processus alatoire, ce qui nous dispense de passer par des calculs fastidieux.

Limite de f (u) quand u 0

1

Nous commenons par rcrire la fonction f 1 :

f1 (u)

= = = E

eiu ln K eiu ln K

(u i) iuF E ST eiu ln ST iuF

ST iu ln ST /K e iuF ST ST sin u ln = E uF K
or lim

sin (u ln ST ) ST K = ln P-ps u0 u K

Donc, par thorme de convergence domine :

30

u0

lim f1 (u)

= E[ =

ST ST ln ] F K

1 E[ST ln ST ] ln K F

1 Pour calculer F E[ST ln ST ], on peut utiliser le thorme de Girsanov pour se ramener aux moments de ln ST calculs auparavant.

En eet :

1 F

E[ST ln ST ] = E[ln ST ] avec E l'esprance sous la probabilit P dnie par : dP dP /Ft = = St rt S e exp
0

0 t

Vs ds + 2

Vs dWs
0

On pose : dWs = dWs de Girsanov).

Vs ds. Sous P, le processus W est un mouvement brownien (Thorme

En remplaant W par W dans l'EDS vrie par (ln St , Vt ), on obtient la dynamique suivante :

d ln St = (r + Vt )dt + Vt dWt 2 dVt = ( + ( )Vt )dt + Vt dZt

Avec Wt et Zt deux mouvements browniens sous la probabilit P tels que < W , Z >t = t En utilisant l'EDS de ln S on obtient : E[ln ST ] = ln F + Pour calculer E[Vt ], nous distinguerons deux cas :
Premier cas : = 0 Dans ce cas, l'EDS de V peut s'crire :

1 2

T 0

E[Vt ]dt.

dVt

( )(

Vt )dt + Vt dZt

Vt dZt

= ( Vt )dt +
avec : = et =
.

Et en utilisant l'esprance de V calcule dans le chapitre 2, on crit :

E[Vt ] = + (V0 )et

En intgrant cette fonction, on arrive enn au rsultat :

f1 (0)

= E[ln ST ] ln K = = ln ln F 1 + K 2
T

E[Vt ]dt
0

F 1 e()T 1 ( )T 1 e()T 1 + + V0 2 K 2 ( ) 2

31

Deuxime cas : = 0 Dans ce cas, l'EDS de V se simplie :

dVt = dt +
et la moyenne de V s'en dduit :

Vt dZt

E[Vt ] = V0 + t
Par suite :

f1 (0)

= ln =

F 1 T E[Vt ]dt + K 2 0 F V0 2 ln + T+ T K 2 4

Remarque : La limite f1 (0) peut exploser quand > 0 et pour des maturits grandes. En plus, un dveloppement limit montre que f1 (u) est quivalente une parabole au voisinage de 0. ie : f1 (u) = f1 (0) + au2 + o(u3 ). La constante a peut prendre des valeurs trs grandes, ce qui donne une parabole trs convexe. Ceci peut poser quelques problmes lors de l'intgration numrique de f1 car on peut avoir f1 (0) de l'ordre de 1030 alors que |f1 (1)| < 1.

En eet, selon le choix des paramtres et des caractristiques de l'option, la forme de la fonction intgrer change. Exemple : Nous avons trac le graphe de f (u) = F f1 (u) K f2 (u) pour deux jeux de paramtres qui montrent cette dicult. La premire forme reste trs lisse, donc facile intgrer, contrairement la deuxime qui a une trs grande pente au voisinage de zro.

Fig.

A.1  Prol de la fonction intgrer pour deux jeux de paramtres.

32

Limite de f (u) quand u 0

2

La deuxime quantit se calcule de faon analogue :

f2 (u)

= = E = E

eiu ln K

(u) iu

eiu ln (ST /K) iu sin u ln ST K u

sin (u ln ST ) ST K = ln P-ps u0 u K Donc, par thorme de convergence domine :

or lim

f2 (0)

= E[ln ST ] ln K ( V0 )(1 eT ) T F + = ln K 2

33

Annexe B
Problmes test

Problme G1 :

13

Minimiser f (x) = 5 Sous les contraintes :

i=1

xi 5
i=1

x2 i
i=5

xi

g1 (x) = 10 2x1 2x2 x10 x11 0 g2 (x) = 10 2x1 2x3 x10 x12 0 g3 (x) = 10 2x2 2x3 x11 x12 0 g4 (x) = 8x1 x10 0 g5 (x) = 8x2 x11 0 g6 (x) = 8x3 x12 0 g7 (x) = 2x4 + x5 x10 0 g8 (x) = 2x6 + x7 x11 0 g9 (x) = 2x8 + x9 x12 0 xi 0, i = 1, . . . , 13 xi 1, i = 1, . . . , 9, 13
Problme G4 :

(B.1)

Minimiser f (x) = 5.3578547x2 + 0.8356891x1 x5 + 37.293239x1 40792.141 3 Sous les contraintes :

avec :

g1 (x) = 92 u(x) 0 g2 (x) = u(x) 0 g3 (x) = 110 v(x) 0 g4 (x) = v(x) 90 0 g5 (x) = 25 w(x) 0 g6 (x) = w(x) 20 0 u(x) = 85.334407 + 0.0056858x2 x5 + 0.0006262x1 x4 0.0022053x3 x5 v(x) = 80.51249 + 0.0071317x2 x5 + 0.0029955x1 x2 + 0.0021813x2 3 w(x) = 9.300961 + 0.0047026x3 x5 + 0.0012547x1 x3 + 0.0019085x3 x4

(B.2)

Problme G6 :

34

Minimiser f (x) = (x1 10)3 + (x2 20)3 Sous les contraintes : (B.3)

g1 (x) = (x1 5)2 + (x2 5)2 100 0 g2 (x) = (x1 6)2 (x2 5)2 + 82.81 0
Problme G8 :

Minimiser f (x) = Sous les contraintes :

sin3 (2x1 ) sin(2x2 ) x3 (x1 + x2 ) 1

(B.4)

g1 (x) = x2 + x2 1 0 1 g2 (x) = x1 1 (x2 4)2 0

Problme G9 :

Minimiser f (x) = (x1 10)2 + 5(x2 12)2 + x4 + 3(x4 11)2 3 +10x6 + 7x2 + x4 4x6 x7 10x6 8x7 5 6 7 Sous les contraintes :

g1 (x) = 2x2 + 3x4 + x3 + 4x2 + 5x5 127 0 1 2 4 g2 (x) = 7x1 + 3x2 + 10x2 + x4 x5 282 0 3 g3 (x) = 23x1 + x2 + 6x2 8x7 196 0 2 6 g4 (x) = 4x2 + x2 3x1 x2 + 2x2 + 5x6 11x7 0 1 2 3
Problme G10 :

(B.5)

Minimiser f (x) = x1 + x2 + x3 Sous les contraintes :

g1 (x) = 1 + 0.0025(x4 + x6 ) 0 g2 (x) = 1 + 0.0025(x4 + x5 + x7 ) 0 g3 (x) = 1 + 0.01(x5 + x8 ) 0 g4 (x) = 100x1 x1 x6 + 833.33252x4 83333.333 0 g5 (x) = x2 x4 x2 x7 1250x4 + 1250x5 0 g6 (x) = x3 x5 x3 x8 2500x5 + 1250000 0

(B.6)

35