Simulation de trajectoires de processus continus

Frédéric Planchet1 et Pierre-E. Thérond2

Résumé. Les processus stochastiques continus sont des outils En effet les données disponibles étant discrètes, l’estimation des
largement employés en finance et en assurance, notamment pour paramètres des modèles en temps continu n’est pas immédiate. De
modéliser taux d’intérêts et cours d’actions. De nombreuses nombreux processus, tels que le modèle de taux proposé par COX,
problématiques amènent à la simulation des trajectoires de tels INGERSOLL et ROSS [1985], n’admettant pas de discrétisation
processus. La mise en œuvre pratique de ces simulations nécessite exacte, l’estimation des paramètres nécessite souvent une
de discrétiser ces processus, d’estimer les paramètres et de générer approximation discrétisée. Les estimations directes à partir du
des nombres aléatoires en vue de générer les trajectoires voulues. modèle discrétisé s’avérant, en général, biaisées lorsque le
Nous abordons successivement ces trois points que nous illustrons processus ne dispose pas d’une version discrète exacte, on est alors
dans quelques situations simples. conduit à utiliser des méthodes indirectes comme la méthode par
inférence indirecte (cf. GIET [2003]) ou encore la méthode des
Mots-clés : Simulation, discrétisation de processus continus, moments efficients pour estimer le modèle.
schémas d’Euler et de Milstein, estimation par inférence indirecte,
générateur du tore mélangé. De plus, de nombreuses problématiques impliquent que l’on soit
capable de simuler l’évolution de cours ou de taux modélisés par
Abstract. Continuous time stochastic processes are useful des processus continus ; la réalisation pratique de telles simulations
models especially for financial and insurance purposes. The nécessitera là encore la discrétisation de ces processus. C’est
numerical simulation of such models is dependant of the time notamment le cas des problématiques de type DFA (Dynamic
discrete discretization, of the parametric estimation and of the Financial Analysis, cf. KAUFMANN, GADMER et KLETT [2001] et
choice of a random number generator. The aim of this paper is to HAMI [2003]) qui consistent en la modélisation de tous les facteurs
provide the tools for the practical implementation of diffusion ayant un impact sur les comptes d’une société d’assurance dans le
processes simulation, particularly for insurance contexts. but d’étudier sa solvabilité ou de déterminer une allocation d’actif
optimale par exemple. Le nombre de variables à modéliser est alors
Keywords: Simulation, Time discrete approximation, Euler and très important : actifs financiers, sinistralité de chaque branche
Milstein Schemes, Indirect inference estimation, mixed torus assurée, inflation, réassurance, concurrence et dépendances entre
random number generator. ces variables. Certaines de ces variables seront modélisées par des
processus à diffusion ; lorsque l’équation différentielle stochastique
1 Introduction (EDS) concernée dispose d’une solution explicite, comme c’est le
cas pour le modèle de VASICEK [1977], la discrétisation s’impose à
La fréquence des variations de cotation des actifs financiers sur les l’utilisateur. Lorsqu’une telle expression n’est pas disponible,
marchés organisés a conduit les financiers à considérer des l’utilisateur pourra se tourner notamment vers les schémas d’Euler
processus stochastiques continus pour modéliser les évolutions de ou de Milstein qui sont des développements d’Itô-Taylor à des
taux d’intérêt comme de cours d’actions. Si les outils ordres plus ou moins importants de l’EDS considérée.
mathématiques sous-jacents peuvent, de prime abord, sembler plus L’approximation sera d’autant plus satisfaisante, au sens du critère
complexes que leurs équivalents discrets, leur usage a notamment de convergence forte (cf. KLOEDEN et PLATEN [1995]), que l’ordre
permis d’obtenir des formules explicites pour l’évaluation d’actifs du développement est élevé.
contingents (cf. les formules fermées de prix d’options européennes
obtenues par BLACK et SCHOLES [1973] lorsque l’actif sous-jacent Par ailleurs, la performance des simulations étant conditionnée par
suit un mouvement brownien géométrique). Ces modèles sont le générateur de nombres aléatoires utilisé, une attention
aujourd’hui utilisés dans de nombreux domaines, notamment en particulière doit être portée à son choix. Nous comparons ici deux
assurance. Leur mise en œuvre pratique nécessite de les discrétiser, générateurs : le générateur pseudo-aléatoire Rnd d’Excel / Visual
que ce soit pour l’estimation des paramètres ou pour la simulation Basic et le générateur quasi-aléatoire obtenu à partir de
des trajectoires. l’algorithme de la translation irrationnelle du tore présenté par
PLANCHET et JACQUEMIN [2003]. L’algorithme du tore s’avère
nettement plus performant que Rnd mais souffre de la dépendance
terme à terme des valeurs générées qui induit un biais dans la
1
Professeur associé de Finance et d’Assurance à l’ISFA et actuaire associé
construction de trajectoires. Pour remédier à ce problème, nous
chez Winter & Associés, fplanchet@[Link]. proposons un générateur du tore mélangé qui conserve les bonnes
2 propriétés de répartition de l’algorithme du tore et peut être utilisé
Chargé de cours en assurance et étudiant en doctorat à l’ISFA et actuaire
pour construire des trajectoires de manière efficace.
chez Winter & Associés, ptherond@[Link].
Institut de Science Financière et d’Assurances (ISFA) – Université
L’objectif du présent article est de proposer au praticien quelques
Lyon 1 - 50, avenue Tony Garnier - 69366 Lyon Cedex 07 - France.
guides méthodologiques lui permettant d’effectuer de manière
Winter & Associés - 9, rue Beaujon - 75008 Paris et 18, avenue Félix efficace des simulations de processus continus, plus
Faure - 69007 Lyon - France.
particulièrement dans le contexte des problématiques d’assurance. où ε est une variable aléatoire de loi normale centrée réduite et δ
Nous abordons donc les trois étapes clé que sont l’estimation des est le pas de discrétisation retenu.
paramètres, la discrétisation du processus et la génération de
nombres aléatoires. 2.2 Discrétisation approximative
2 Discrétisation de processus continus Lorsque la discrétisation exacte n’existe pas, il convient de se
tourner vers des approximations discrètes du processus continu
Prenons le cas d’un processus défini par l’EDS suivante : sous-jacent. Les schémas d’Euler et de Milstein sont les procédés
 dX t = µ ( X t , t ) dt + σ ( X t , t ) dBt
de discrétisation les plus répandus. Tous deux sont des
 (1) développements d’Itô-Taylor de l’équation (2) à des ordres
 X 0 = x différents (ordre 1 pour Euler, ordre 2 pour Milstein).
où B est un mouvement brownien standard.
Dans la suite, nous ferons référence au critère de convergence forte
La mise en oeuvre pratique de ce processus va nécessiter la pour classer les procédés de discrétisation. Une discrétisation
~
discrétisation de celui-ci. Pour cela on met l’équation (1) sous sa approximative X converge fortement vers le processus continu X
forme intégrale : lorsque l’erreur commise sur la valeur finale (à la date T) de la
t t trajectoire obtenue par le processus discrétisé est en moyenne
X t = x + µ ( X s, s ) ds + σ ( X s, s ) dBs .
∫0 ∫0 (2) asymptotiquement négligeable, i. e. lorsque
~δ
∀ T > 0, lim E  X T − X T  = 0 (6)
Si le processus considéré ne dispose pas d’une discrétisation δ → 0  
exacte, un développement d’Itô-Taylor de l’équation (2) nous La vitesse de convergence de l’équation (6) nous permet
permet de disposer d’une version discrétisée approximative (voir le d’introduire un ordre entre les procédés de discrétisation. Ainsi le
paragraphe 2.2). Cette approximation est d’autant plus précise que ~
processus discrétisé X converge fortement à l’ordre γ vers le
le développement intervient à un ordre élevé. processus X si :
∃ K > 0, ∃δ0 > 0, ∀δ ∈ ]0, δ0 ] E  X T − X T  ≤ K δ .
~δ γ
2.1 Discrétisation exacte   (7)

La simulation d’un processus d’Itô pourra être effectuée On notera que ce critère ne fait pas référence à la convergence
directement (sans erreur de discrétisation) dès lors que celui-ci uniforme.
admet une discrétisation exacte. Rappelons la définition d’une
discrétisation exacte. 2.2.1 Schéma d’Euler
Le procédé de discrétisation d’Euler consiste en l’approximation
Définition : Un processus (X~kδ )k∈[1;T / δ ]
est une discrétisation ~
du processus continu X par le processus discret X défini, avec les
~
exacte du processus X si ∀ δ > 0, ∀ k ∈ [1; T / δ ], X k δ =loi X k δ .
mêmes notations que précédemment, par :
~ ~ ~
( ) ~
( )
X t + δ = X t + µ X t, t δ + σ X t, t δ ε . (8)
Un processus X admet une discrétisation exacte dès lors que l’on
KLOEDEN et PLATEN [1995] prouvent que sous certaines conditions
peut résoudre explicitement l’EDS qui lui est associée. C’est
de régularité, le schéma d’Euler présente un ordre de convergence
notamment le cas du mouvement brownien géométrique retenu par
forte de 0,5.
BLACK et SCHOLES [1973] pour modéliser le cours d’une action ou
encore celui du processus d’Ornstein-Ulhenbeck retenu par
Par exemple dans le modèle de Cox, Ingersoll et Ross (CIR), le
VASICEK [1977] pour modéliser le taux d’intérêt instantané r :
taux d’intérêt instantané est solution de l’EDS :
drt = a (b − rt ) dt + σ dBt .
drt = a (b − rt ) dt + σ rt dBt .
(3)
(9)
On rappelle que l’équation (3) exprimée dans l’univers historique
Aussi le processus discret ~
r déterminé par le schéma d’Euler peut
ne suffit pas à définir le prix des zéro-coupons, et doit être
s’écrire :
complétée par la forme de la prime de risque. Celle-ci est choisie
constante dans le modèle de VASICEK, ce qui permet de conserver rt + δ = r~t + a ( b − r~t ) δ + σ ~
~ rt * δ ε . (10)
la même forme de la dynamique dans l’univers corrigé du risque.
L’équation (2) s’écrit dans ici :
2.2.2 Schéma de Milstein
( )
t
rt = r0 e− at + b 1 − e− at + σ e− at eas dBs . ∫0 (4) Le schéma de Milstein est obtenu en allant plus avant dans le
~
développement d’Itô-Taylor. Le processus discret X est alors
Les propriétés de l’intégrale d’une fonction déterministe par
défini par :
rapport à un mouvement brownien conduisent à la discrétisation
exacte : ~ ~ ~
( ) ~
( )
X t +δ = X t + µ X t , t δ + σ X t , t δ ε

( ) ( ) ( ) ( )
− 2aδ ~ ~ (11)
− aδ − aδ 1− e (5) σ X , t σ Xt, t
rt + δ = rt e + b 1− e +σ ε, + x t δ ε2 − 1
2a 2
 ( )
~
où σ x X t , t désigne la dérivée par rapport au premier argument de peuvent faire intervenir plus d’une variable aléatoire ce qui signifie

( )
~ des temps de simulation plus importants.
la fonction σ(.,.) évaluée en X t , t .
En outre, le schéma de Milstein pose des problèmes pratiques de
Ce procédé de discrétisation présente, en général, un ordre de mise en œuvre dès lors que l’on s’intéresse à des phénomènes de
convergence forte de 1. Remarquons que si la volatilité σ est une dimension supérieure ou égale à 2, i. e. dès que l’on cherche à
constante, les procédés de discrétisation d’Euler et de Milstein simuler des vecteurs aléatoires.
conduiront à la même discrétisation, c’est le cas pour les modèles
de Vasicek et de Hull et White (cf. HULL [1999]). En revanche Le graphique suivant permet de comparer les évolutions moyennes
pour le modèle de CIR, il vient : de la diffusion définie par (3) selon le schéma de discrétisation

)
retenu pour les paramètres :
(
2
σ
~ rt + a (b − ~
rt + δ = ~ rt ) δ + σ ~
rt * δ ε + δ ε2 − 1 . (12) r0 = 4% b = 5% a = 0,5 σ = 10%
4
Ce graphique a été obtenu en faisant la moyenne de 10 000
En développant à des degrés supérieurs l’équation (2), il est
simulations du taux d’intérêt générées.
possible d’obtenir des processus discrétisés d’ordre de convergence
plus élevé. Toutefois ils nécessiteront des calculs plus nombreux et

5,25%

5,00%

4,75%

4,50%

Schéma d'Euler

4,25% Discrétisation exacte

Espérance

4,00%
0 5 10 15 20 25 30 35 4

Fig. 1. Evolution moyenne du taux modélisé par Vasicek selon le procédé de discrétisation retenu

Les discrétisations exacte et selon le schéma d’Euler sont δr

dynamique BCt + δ = BCt e t . A partir des taux simulés pour
relativement proches graphiquement. L’écart est en particulier
négligeable lorsqu’on s’intéresse aux évolutions moyennes d’un effectuer le graphique 1, on obtient le graphique suivant.
bon de capitalisation de prix BCt à l’instant t, qui évolue selon la

5,5

5,0

4,5
Schéma d'Euler
4,0
Discrétisation exacte
3,5
Espérance
3,0

2,5

2,0

1,5

1,0
0 5 10 15 20 25

Fig. 2. Evolution moyenne d’un bon de capitalisation dont le taux évolue

selon le modèle de Vasicek selon le procédé de discrétisation retenu
 approximative de (1), le schéma d’Euler est souvent utilisé pour
3 Estimation des paramètres servir de modèle auxiliaire.

L’estimation des paramètres est une étape délicate dans la Cette méthode conduit à des estimations bien plus précises que
simulation de trajectoires d’un processus continu car elle peut être celles obtenues à partir d’estimations « naïves » des paramètres du
l’origine d’un biais. En effet, le praticien peut se voir confronté à modèle auxiliaire sur les données observées.
deux problèmes :
le processus n’admet pas forcément de discrétisation GIET [2003] étudie l’impact du choix du procédé de discrétisation
exacte, utilisé pour fournir le modèle auxiliaire sur l’estimation par
la variable modélisée n’est pas toujours directement inférence indirecte et montre qu’utiliser un procédé d’ordre de
observable. convergence plus élevé (comme la schéma de Milstein par rapport
au schéma d’Euler) permet de réduire considérablement le biais sur
En effet si le processus considéré n’admet pas de discrétisation l’estimation des paramètres.
exacte, il sera impossible d’estimer les paramètres du modèle par la
méthode du maximum de vraisemblance. Par ailleurs il est souvent Toutefois, que ce soit dans le cas d’une estimation par maximum
impossible d’exprimer de manière exacte les densités transitoires de vraisemblance ou par inférence indirecte, on doit observer que
entre deux observations. Il faudra alors se tourner vers des dans le cas d’un modèle de taux tel que celui de Vasicek,
méthodes simulées telles que l’inférence indirecte pour estimer les l’estimation des paramètres de la dynamique d’évolution du taux
paramètres. court dans l’univers historique ne suffit pas et doit être complétée
par l’estimation de la prime de risque, qui constitue un point délicat
3.1 Le biais associé à la procédure d’estimation (cf. LAMBERTON et LAPEYRE [1997] sur ce sujet).

On reprend ici la cas simple de la modélisation du taux à court 3.3 Les méthodes d’estimation ad hoc
terme par un processus d’Ornstein-Ulhenbeck, comme dans
l’équation (3) ; la discrétisation exacte (5) conduit à un processus En pratique, dans les problèmes assurantiels, la grandeur modélisée
auto-régressif d’ordre 1, AR(1). par un processus continu n’est, le plus souvent, pas la grandeur
d’intérêt : typiquement, on modélise le taux court, par exemple par
L’estimation des paramètres du modèle s’effectue classiquement en un processus de diffusion, mais dans le cadre d’une problématique
régressant la série des taux courts sur la série décalée d’une d’allocation d’actif pour un régime de rentiers, il sera déterminant
période, que l’on écrit sous la forme usuelle : que le modèle représente correctement les prix des obligations
d’échéances longues. De manière équivalente, on est amené à
Y = α + β X + σ1 ε , (13) estimer les paramètres du modèle pour représenter correctement la
avec courbes des taux zéro-coupon, qui est dans une relation bijective
avec la courbe du prix des obligations zéro-coupon :
α 2 ln β
a = − ln β , b = et σ 2 = σ12 2 .
P(0, T ) = exp {− T R (0, T )} .
(14)
1−β β −1 (16)
L’estimateur du maximum de vraisemblance (EMV) coïncide avec L’idée est alors d’estimer les paramètres pour minimiser une
l’estimateur des moindres carrés ; on obtient ainsi les estimateurs distance (par exemple la distance quadratique) entre les prix prédits
suivants des paramètres du modèle d’origine : par le modèle et les prix observés sur le marché. Cette méthode est
αˆ 2 ln βˆ notamment utilisée dans FARGEON et NISSAN [2003].
aˆexact = − ln β̂ , bˆexact =
2
et σˆ exact = σˆ 2 2 . (15)
1− βˆ ˆ
β −1 De plus, en pratique on pourra s’interroger sur les paramètres que
l’on estime et les paramètres que l’on fixe arbitrairement compte
On constate en particulier que ces estimateurs, s’ils sont EMV, sont tenu de la connaissance que l’on peut avoir du contexte par
biaisés. Ceci peut être gênant, notamment dans le cas où les ailleurs. Ainsi, dans l’exemple évoqué, l’estimation simultanée des
paramètres du modèle ont une interprétation dans le modèle, paramètres a, b et σ conduit à une courbe des taux quasi
comme par exemple b dans le cas présent, valeur limite du taux déterministe (σ est petit), ce qui peut apparaître irréaliste. On est
court. On voit sur cet exemple qu’on sera donc conduit à mettre en ainsi conduit à fixer arbitrairement σ, puis à estimer les deux
œuvre des procédés de réduction du biais. paramètres restants.

3.2 Principe de l’estimation par inférence indirecte Dans cette approche il convient également d’être attentif au choix
effectué quand à la paramétrisation du modèle : courbe des taux
Cette méthode est utilisée lorsque le processus n’admet pas de zéro-coupon ou courbe des prix des zéro-coupon ; en effet, la
discrétisation exacte ou que sa vraisemblance, trop complexe, ne fonction de correspondance entre ces deux grandeurs accentue les
permet pas d’implémenter la méthode du maximum de écarts de courbure et le choix de l’une ou l’autre des
vraisemblance. Elle consiste à choisir le paramètre θ̂ qui minimise paramétrisations peut conduire à des résultats sensiblement
(dans une métrique à définir) la distance entre l’estimation d’un différents. Ce point est développé sur le plan théorique par
modèle auxiliaire sur les observations et l’estimation de ce même RONCALLI [1998] et est illustré sur un exemple dans la section
modèle auxiliaire sur les données simulées à partir du modèle de suivante.
base pour θ = θˆ . Le modèle auxiliaire étant une discrétisation
Au surplus, on notera que l’exploitation directe de données telles
que le prix des zéro-coupon (ou les taux zéro-coupon) évite
l’estimation de la prime de risque ; celle-ci est en effet incluse dans et d’autre part à partir des prix des obligations qui s’en déduisent.
ces prix de marché et conduit ainsi à travailler naturellement dans Cette courbe des taux est construite à partir des prix de marché et
l’univers corrigé du risque. L’estimation directe des paramètres sur elle inclut les primes de risque.
des historiques de taux court nécessite, on l’a vu, l’estimation
(délicate) de la prime de risque. Ce point est illustré par exemple Estimation ad hoc sur le prix des zéro-coupons : cette procédure
dans LAMBERTON et LAPEYRE [1997]. consiste à déterminer par une méthode de type « moindres carrés »
les paramètres du modèle de Vasicek qui permettent de représenter
Dans les illustration présentées ci-après, l’exploitation de données au mieux les prix des zéro-coupons. Cette méthode a été mise en
incluant les primes de risque est privilégiée. ~
œuvre en estimant les trois paramètres a, b et σ puis en fixant σ
et en n’estimant plus que le deux autres paramètres.
3.4 Illustration dans le cas du modèle de Vasicek
Comme on peut le voir dans le tableau 1, l’estimation des
A partir de la courbe des taux publiée par l’Institut des Actuaires,
paramètres du modèle de Vasicek donne des résultats très
nous avons estimé les paramètres du modèle de Vasicek selon
différents selon l’approche retenue :
l’approche ad hoc décrite supra. L’estimation est menée d’une part
directement à partir des valeurs des taux zéro-coupon de la courbe

EMC sur le prix des zéro-coupon EMC sur les taux zéro-coupon
Paramètres
a, b et σ a et b a et b a, b et σ a et b a et b
estimés
a= 0,2293 0,2353 0,5032 0,3951 0,3253 0,4232

b= 0,0572 0,0760 0,0701 0,0624 0,0745 0,0869

r0 = 0,0207 0,0207 0,0207 0,0207 0,0207 0,0207

σ= 0,0000 0,0500 0,1000 0,0000 0,0500 0,1000

Tab. 1. Paramètres du modèle de Vasicek estimés selon les techniques d’estimation des paramètres

Le prix P(t) en 0 d’un zéro-coupon d’échéance t est donné par la La figure 3, qui indique le prix des zéro-coupons, en fonction de
formule classique : leur échéance, selon les différentes méthodes d’estimation des
  paramètres, nous permet d’observer que l’estimation des
1   ~ σ 2 
( )
2
 ~ σ paramètres par le critère des moindres carrés sur les taux zéro-
P(t) = exp  − t *  b − −  b − − r0  1 − e − a t
  2a at   2a  coupon s’avère moins performant que la méthode ad hoc si la
   
(17) variable d’intérêt est le prix des obligations zéro-coupon.
 
−
σ2
4a
(
1− e −a t 2   
 
)
  
1,0

Observation

0,8 Estimation MCO sur les prix

Estimation MCO sur les taux

0,6

0,4

0,2

0,0
0 5 10 15 20 25 30

Fig. 3. Prix des zéro-coupons en fonction de leur échéance

et de la méthode d’estimation des paramètres
En revanche, cette technique d’estimation des paramètres permet espérance de la courbe des taux originelle que ceux obtenus par la
de générer des taux zéro-coupon nettement plus proches en méthode ad hoc :

0,07

0,06

0,05

0,04

0,03

Observation
0,02
Estimation MCO sur les taux

0,01 Estimation MCO sur les prix

0,00
0 5 10 15 20 25 30

Fig. 4. Taux instantanés espérés en fonction de leur échéance

et de la technique d’estimation des paramètres

Si les taux obtenus par la méthode ad hoc sur le prix des zéro- l’algorithme de Box-Muller permet à partir de deux variables
coupons sont, en espérance, éloignés des taux observés, ils uniformes indépendantes sur [0,1] de générer deux variables
permettent néanmoins d’approcher le prix des zéro-coupons avec indépendantes N(0,1). En effet, si U1 et U2 sont des v.a. U[0,1]
une bonne précision. indépendantes alors en posant :
 X = − 2 ln U cos (2πU )
Ainsi, en pratique, on sera conduit à privilégier pour une  1 1 2
 (18)
problématique d’allocation d’actifs d’un régime de rentes  X 2 = − 2 ln U1 sin (2πU 2 )
l’approche ad hoc sur le prix des zéro-coupons.
Les variables aléatoires X1 et X2 sont indépendantes et suivent
une loi gaussienne centrée réduite. Cette technique (exacte) est
4 Génération des trajectoires toutefois relativement longue à mettre en œuvre et nécessite
La génération des trajectoires passe nécessairement par la l’indépendance des réalisations uniformes générées. AUGROS et
génération de nombres aléatoires. De manière pratique, il s’agit MORÉNO [2002] présentent divers algorithmes pour approximer
de générer des réalisations de variable aléatoire de loi uniforme l’inverse de la fonction de répartition de la loi normale centrée
sur le segment [0, 1]. En effet, si u est une telle réalisation, réduite. Dans la suite, nous utiliserons l’algorithme de Moro (cf.
PLANCHET et JACQUEMIN [2003]) qui allie rapidité et précision.
F −1(u) peut s’apparenter à une réalisation d’une variable
aléatoire de fonction de répartition F. La technique d’inversion 4.1 Deux générateurs de nombres aléatoires : Rnd
de la fonction de répartition permet ainsi à partir de réalisations
de variables uniformes, d’obtenir de simuler des réalisations et l’algorithme du tore
d’autres variables aléatoires. Lorsqu’on ne dispose pas de
formule explicite pour F-1, on utilisera des algorithmes 4.1.1 Le générateur implémenté dans Excel / Visual Basic
d’approximation de cette fonction ou des algorithmes spécifiques
à la loi que l’on souhaite traiter Par exemple, pour simuler une Le générateur implémenté dans Excel (Rnd) est un générateur
réalisation d’une loi de Poisson de paramètre λ, on générera une congruentiel, c’est à dire un générateur périodique issu d’une
suite (Vi ) de réalisations de variable aléatoire exponentielle de
valeur initiale (on parle également de « graine » du générateur).
Changer de valeur initiale permet de changer de suite de
 n +1  nombres. Le lecteur se réfèrera à PLANCHET et JACQUEMIN [2003]
paramètre 1. En effet N = Inf  n

∑ Vi > λ  suit une loi de pour plus d’informations sur l’implémentation informatique de
i =1  tels générateurs.
Poisson de paramètre λ.
4.1.2 La translation irrationnelle du tore
Les modélisations retenues en finance et en assurance faisant Ce générateur multidimensionnel donne à la n-ème réalisation de
souvent intervenir des mouvements browniens, il est nécessaire la d-ème variable aléatoire uniforme à simuler la valeur un :
de simuler des réalisations de variables aléatoires N(0,1). On ne
dispose pas de formule exacte de l’inverse de la fonction de
répartition inverse de la loi normale centrée réduite, mais
[
u n = n pd − n pd , ] (19)
où pd est le d-ème nombre premier et [.] désigne l’opérateur annexe permettent d’apprécier la supériorité de répartition de la
partie entière. suite générée par le tore par rapport à celle générée par Rnd.

4.1.3 Quelques éléments de comparaison

Les tests statistiques élémentaires (test d’adéquation du χ2, de
Kolmogorov-Smirnov et d’Anderson-Darling) rappelés en

Fig. 5. Répartition et tests statistiques du générateur Rnd

Fig. 6. Répartition et tests statistiques du générateur du tore*

Par ailleurs, pour comparer les performances des deux générateurs,  

 σ2  
nous avons étudié la précision de l’estimation empirique du prix St + δ = St exp   r − δ+σ δε , (21)
d’un call européen sur une action modélisée par un mouvement  2 
   
brownien géométrique sous la probabilité risque-neutre :
où ε est une variable aléatoire de loi N(0,1).
~
dSt = r St dt + σ St dBt , (20)
~ Par ailleurs, nous disposons d’une formule fermée pour le prix en t
où r est le taux sans risque et B un mouvement brownien sous la d’un Call européen de prix d’exercice K et d’échéance T sur ce
mesure risque-neutre. Rappelons que, sous les hypothèses du titre :
modèle de Black et Scholes, la probabilité risque-neutre est
l’unique probabilité sous lesquels les prix actualisés sont des Ct (S, K, T) = St N(d1 ) − K e−r(T − t) N(d 2) , (22)
martingales. La densité de Radon-Nikodym de cette probabilité par où :
rapport à la probabilité historique est donnée par le théorème de
Girsanov (cf. AUGROS et MORÉNO [2002]). St  σ 2 
ln + r+ (T − t )
K  2 
Ce processus dispose d’une discrétisation exacte : d1 =   ,
σ T −t
 (T − t )
d 2 = d1 − σ T − t ,
N désigne la fonction de répartition de la loi normale
Cˆt ( S, K, T, n ) =
e −r
n ∑
n
[sTi − K ]
+
. (23)
i =1
centrée réduite.
L’erreur relative d’estimation peut donc s’écrire :
Nous pouvons donc mesurer la performance des générateurs par le
biais de l’erreur d’estimation du prix de l’option considérée. En Cˆ ( S, K, T, n ) − Ct ( S, K, T )
ρ= t . (24)
effet, si sTi est la prix du titre considéré à la date T dans la i-ème Ct ( S, K, T )
simulation, l’estimateur naturel du prix de l’option en t est : Le graphique suivant présente l’évolution de ρ selon le nombre de
simulations effectuées pour une action de volatilité 20 % et un Call
d’échéance 6 mois.
1%

0%
0 10 000 20 000 30 000 40 000 50 000 60 000 70 000 80 000 90 000 100 000

-1%

Algorithme du tore (p = 5)

Générateur Rnd

-2%

-3%

Fig. 7. Erreur d’estimation en fonction du nombre de simulations

L’estimation à partir des valeurs générées par Rnd est biaisée de volatilité 20 % et avons comparé l’évolution moyenne du cours
manière systématique de prés de 0,2 % alors celle effectuée à l’aide estimée à partir des simulations et l’évolution moyenne théorique.
du tore converge rapidement vers la valeur théorique : à partir de
7 000 simulations, l’erreur relative est inférieure à 0,1 %. En effet,si sti est le cours du titre à la date t dans la i-ème

Il ressort des ces différentes comparaisons que le générateur du

simulation, l’estimateur empirique du cours moyen à la date t, st
tore semble plus performant que Rnd. Son utilisation connaît est donné par :
toutefois un certain nombre de limites. 1 n i
st = ∑ s .
n i =1 t
(25)
4.2 Limites de l’algorithme du tore Le graphique suivant présente les deux trajectoires obtenues.
Les valeurs générées par le tore ne sont pas indépendantes terme à
terme, ceci peut générer des erreurs non négligeables. Par exemple
nous avons généré 10 000 trajectoires, avec un pas de discrétisation
de 1, sur 20 ans d’un mouvement brownien géométrique de
 275

250
Espérance
Moyenne empirique à partir des trajectoires générées par le tore (p=5)
225

200

175

150

125

100
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Fig. 8. Trajectoires espérée et moyenne des 10 000 trajectoires simulées

par le tore d’un mouvement brownien géométrique

A partir de t = 2, les trajectoires générées à partir de l’algorithme trajectoires. Ce biais s’explique par la dépendance terme à terme
du tore sont en moyenne très en dessous de la trajectoire moyenne. des valeurs générées par le tore que nous permet d’observer le
Cet algorithme n’est donc pas utilisable en l’état pour générer des corrélogramme de la suite générée par ce générateur.

100%

75%

50%

25%

0%
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

-25%

-50%

Fig. 9. Corrélogramme de la suite générée par le tore (p = 5)

Rappelons que le h-ème terme du corrélogramme ρh s’écrit : test consiste à prendre des listes de quatre chiffres tirés
aléatoirement, de manière à observer les « combinaisons de
n valeurs » puis à effectuer un test d’adéquation du χ2 pour voir si les
∑ (uk − u )(uk + h − u )
k =1
fréquences observées correspondent aux fréquences théoriques. On
ρh = , (26) distingue ainsi cinq cas :
n
∑ (uk − u ) 2 les quatre chiffres sont tous différents,
la liste contient une et une seule paire,
k =1
la liste contient deux paires,
où u désigne la moyenne empirique de la suite u. la liste contient un brelan (trois chiffres identiques),
la liste contient un carré (quatre chiffres identiques).
Le « test du poker » permet également de mettre en évidence cette
faille de l’algorithme du tore. L’idée de ce test est de comparer les En constituant les listes en prenant les chiffres dans l’ordre où ils
fréquences théoriques des mains au poker avec les fréquences sont simulés, le tableau 2 résume les fréquences suivantes.
observées sur les simulations effectuées. De manière pratique, ce
 Tab. 2. Test du poker sur 4 000 réalisations générées par l’algorithme du tore

Les p-valeurs des tests d’adéquation du χ2 sont proches de 1,

um = uφ (n) , (27)
toutefois quel que soit p, le tore ne conduit jamais à l’obtention
d’un carré, ni d’un brelan, ni d’une double paire alors que ces trois où :
situations représentent 6,4 % des cas.
φ(n) = [α * N * u~ + 1] , (28)
4.3 Générateur du tore mélangé où :
Pour contourner le problème posé par la dépendance terme à terme [.] désigne l’opérateur partie entière,
des valeurs générées par l’algorithme du tore, nous proposons une α ≥ 1,
adaptation qui consiste à « mélanger » ces valeurs avant de les u~ est la réalisation d’une variable aléatoire de loi
utiliser. uniforme.

4.3.1 Descriptif de l’algorithme Le facteur α de l’équation (28) a pour vocation de réduire le

nombre de tirages qui donneraient lieu au même indice et donc au
Notons (un ) la suite générée par le nombre premier p. Au lieu même nombre aléatoire. En effet, plus α est grand plus la
d’utiliser la nombre un lors du n-ème tirage de la loi uniforme sur probabilité de tirer deux fois le même nombre aléatoire est faible.
[0, 1], nous proposons d’utiliser um où m est choisi de manière Dans la pratique, α = 10 est satisfaisant.
aléatoire dans N.
4.3.2 Choix de la procédure de mélange
Le générateur ainsi obtenu présente les mêmes bonnes
Pour la génération de u~ , nous avons retenu le générateur Rnd ou
caractéristiques globales que l’algorithme du tore sans la
dépendance terme à terme. Il nécessite toutefois davantage de tout autre générateur congruentiel ayant une période importante. En
temps de simulation du fait du tirage de l’indice m. Nous effet, utiliser l’algorithme du tore pour effectuer le mélange ne
proposons l’algorithme suivant lorsque l’on souhaite générer N réduit pas la corrélation observée comme le montre le
réalisation de variables uniformes : corrélogramme suivant.

100%

75%

50%

25%

0%
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

-25%

-50%

Fig. 10. Corrélogramme de la suite générée par le tore (p1 = 5) mélangé par le tore (p2 = 19)
 100%

75%

50%

25%

0%
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

-25%

-50%

Fig. 11. Corrélogramme de la suite générée par le tore (p1 = 5) mélangé par Rnd

En revanche l’utilisation de Rnd permet de réduire De plus, comme le montre le tableau 3, le tore mélangé par Rnd
considérablement les corrélations. La figure 11 nous permet de satisfait de manière satisfaisante le test du poker puisque toutes les
constater que le mélange a pratiquement fait disparaître la « mains » sont représentées dans des proportions proches des
corrélation terme à terme, le générateur du tore mélangé fréquences théoriques.
conservant néanmoins les propriétés de bonne répartition globale et
de rapidité de convergence du tore.

Tab. 3. Test du poker sur 40 000 réalisations générées par le tore mélangé par Rnd

275

250
Espérance
Moyenne empirique à partir des trajectoires générées par le tore mélangé (p=5)
225

200

175

150

125

100
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Fig. 12. Trajectoires espérée et moyenne des 5 000 trajectoires simulées

par le tore mélangé d’un mouvement brownien géométrique
La figure 12 nous permet de vérifier que les trajectoires générées à [15] Vasicek O. [1977] « An equilibrium characterization f the term
partir du tore mélangé sont satisfaisantes. structure ». Journal of financial Economics, vol. 5, 177-188.

5 Conclusion Annexe : Tests d’adéquation à une loi

Dans cet article nous avons présenté les trois étapes indispensables Le lecteur trouvera un descriptif plus détaillé de ces tests dans
à la génération pratique de trajectoires de variables modélisées par SAPORTA [1990] pour les deux premiers et dans PARTRAT et
des processus continus : l’estimation des paramètres, la BESSON [2004] pour le dernier.
discrétisation du processus et la génération de nombres aléatoires.
Test du χ2
L’étape de discrétisation conduit à arbitrer entre précision et temps
de calcul à moins qu’une discrétisation exacte peu coûteuse en Soit X une variable aléatoire à valeurs dans E de loi PX inconnue et
calculs soit disponible. Lorsque ce n’est pas le cas, on préfèrera le P0 une loi connue sur E. Soit A1, …,Ad une partition de E telle que
schéma de Milstein au schéma d’Euler car il demande le même π k = P0 [ Ak ] > 0 pour tout k. Nous disposons d’un n-échantillon
nombre de tirages de nombres aléatoires pour une meilleure indépendant et identiquement distribué (X1, …,Xn) de X. Soit Nk le
précision. nombre de variables aléatoires Xi dans Ak. Si PX = P0 (H0)
considérons la statistique D2 définie comme suit :
L’estimation des paramètres doit faire l’objet d’une attention d ( Nk − n πk )2 .
particulière du fait des biais auxquels conduirait une estimation
« naïve ». Il convient de préférer la méthode d’estimation par
2
D = ∑
k =1 n πk
(29)

inférence indirecte à une estimation directe souvent porteuse de

biais. Le principe de l’estimation ad hoc permettra également
d’estimer des paramètres satisfaisants lorsque la variable d’intérêt D2 est asymptotiquement distribuée comme une variable de χ 2d −1 .
n’est pas directement la variable simulée, ce qui est le plus souvent La p-valeur de ce test est donc donnée par :
le cas dans les problématiques d’assurance.
[2 2
αˆ = P χ d −1 > D . ] (30)
Enfin nous recommandons l’utilisation du tore mélangé pour toute
construction pas-à-pas de trajectoires. Ce générateur, simple à
mettre en place, donne des résultats très satisfaisants, notamment Test de Kolmogorov-Smirnov
comparé au Rnd d’Excel.
Si Fn* est la fonction de répartition empirique d’un n-échantillon
REFERENCES d’une variable aléatoire de fonction de répartition F, alors la
[1] Augros J.C., Moreno M. [2002] Les dérivés financiers et *
statistique Dn = Sup Fn (x) − F(x) est asymptotiquement
d’assurance. Paris : Economica.
[2] Black F., Scholes M. [1973] « The pricing of options and corporate distribuée comme suit :

] ∑ (−1) exp {− 2 k y }.
liabilities ». Journal of Political Economy, vol. 81, n°3, 637-654.
[3] Cox J.C., Ingersoll J.E., Ross S.A. [1985] « A theory of the term
structure of interest rates ». Econometrica, vol. 53, 385-407.
P [ nDn <y →
∞
k 2 2
(31)
−∞
[4] Fargeon L., Nissan K. [2003] Recherche d’un modèle actuariel
d’analyse dynamique de la solvabilité d’un portefeuille de rentes
viagères. Mémoire d’actuariat, ENSAE. La statistique Dn est indépendante de F et le résultat (31) nous
[5] Giet L. [2003] « Estimation par inférence indirecte des équations de permet de tester :
diffusion : l’impact du choix du procédé de discrétisation ».  H 0 : F(x) = F0(x)
Document de travail n°03A15 du GREQAM. 
[6] Hami S. [2003] Les modèles DFA : présentation, utilité et
 H1 : F(x) ≠ F0(x) (32)
application. Mémoire d’actuariat, ISFA. 
[7] Hull J. [1999] Options, futures and other derivatives, 4th edition.  α
Prentice-Hall.
[8] Kaufmann R., Gadmer A., Klett R. [2001] « Introduction to Dynamic
Financial Analysis ». ASTIN Bulletin, vol. 31, n°1, 213-249. Test d’Anderson-Darling
[9] Kloeden P., Platen E. [1995] Numerical Solution of Stochastic
Differential Equations, 2nd Edition. Springer-Verlag.
[10] Lamberton D., Lapeyre B. [1997] Introduction au calcul
Ce test repose sur l’écart d’Anderson-Darling An2 défini comme
stochastique appliqué à la finance, 2e édition, Paris : Ellipses. suit :
[11] Partrat Ch., Besson J.L. [2004] Assurance non-vie. Modélisation,
simulation. Paris : Economica.
∞
(F (x) − F(x) ) dF ( x ) .
* 2

[12] Planchet F., Jacquemin J. [2003] « L’utilisation de méthodes de An2 = n ∫0 F(x) (1 − F(x) )
n
0
(33)
simulation en assurances ». Bulletin Français d’Actuariat, vol. 6,
n°11, 3-69.
[13] Roncalli T. [1998] La structure par terme des taux zéro :
modélisation et implémentation numérique. Thèse de l’université Cet écart permet de tester (32), l’expression opérationnelle de An2
Montesquieu-Bordeaux IV. étant :
[14] Saporta G. [1990] Probabilités, analyse des données et statistique.
Paris : Editions Technip. An2 = −n −
1 n
∑
n i =1
{ [
(2i − 1) ln F(x(i)) + ln 1 − F(x(n − i +1)) (34) ]}
où x(i) désigne la i-ème plus petite réalisation de X dans
est supérieur à une valeur que la variable aléatoire An2 a une
l’échantillon.
probabilité α de dépasser.
On rejettera H0 si
ACKNOWLEDGEMENTS
−n−
1 n
∑
n i =1
{ [
(2i − 1) ln F0(x(i)) + ln 1 − F0(x(n − i +1)) ]} (35)
Les auteurs tiennent à remercier Pierre Devolder pour ses
remarques qui nous ont aidées dans la finalisation de cet article.