Hasbnclic 708
Hasbnclic 708
Florian Ielpo1
24 février 2008
1
Dexia Group, 7/11 Quai André Citroën, 75015 Paris, Centre d’Economie de la Sorbonne
- Antenne de Cachan, Avenue du Président Wilson, 94230 Cachan. E-mail : [Link]@clf-
[Link]
2
Table des matières
3
4 TABLE DES MATIÈRES
Bibliographie 171
Introductions
”S o if you do not accept the Gaussian distribution (i.e. if you have some ethics) AND do not
”value” options in a axiomatized top-down fashion (i.e. only price them as some informed and
temporary guess), then YOU ARE NOT USING THE BLACK SCHOLES FORMULA, but one
of the modifications of the one started by Bachelier (the latest contributer being Ed Thorp’s).
They did not ground their formula in the Gaussian.”
Nassim Nicholas Taleb1
Il s’agit dans un premier temps de revenir sur le modèle linéaire gaussien, dans sa ver-
sion univariée et multivariée. On présentera quelques questions simples liées à l’inférence
statistique de ces modèles ainsi qu’à l’usage qui peut en être fait : expliquer la dynam-
qiue des séries économiques/financières et permettre la mise en oeuvre de prévisions
1
[Link]
7
8 TABLE DES MATIÈRES
Il s’agit ensuite de présenter la base de la théorie des séries temporelles : modèle ARMA,
GARCH et modèles à facteurs. Là encore, la principale motivation sera l’inférence ef-
ficace ainsi que la prévision.
Ces notes de cours s’appuient sur un certain nombre d’ouvrages de statistiques bien
connus ainsi que sur d’autres notes de cours, qui seront citées à chaque fois. Nous y
renvoyons un lecteur soucieux de dépasser le niveau de cette introduction. La partie
consacrée au modèle linéaire gaussien est grandement inspirée de Greene (2002). La par-
tie consacrée à l’étude des séries temporelles est principalement inspirée de Cochrane
(2005).
Rappels de mathématiques et
probabilité
Cette première partie a pour but de revenir sur un certain nombre de concepts et
techniques nécessaires pour comprendre et implémenter les différentes méthodes de
l’économétrie de la finance. Il s’agit principalement de revenir sur un certain nombre
de concepts de probabilités dans un premier temps (définition d’une variable aléatoire,
de ses moments et des distributions qu’il est possible de lui affecter). Il sera ensuite
question de revenir sur les concepts de convergence (presque sure, en probabilité et
en loi), afin d’introduire la Loi des Grands Nombres (LGN hereafter) et le Théorème
Central Limite (TCL). Enfin, on finira par quelques éléments de calculs matriciel.
Définition 1.1.1 (Notion de tribu). On dit qu’une partie A de P(Ω) est une tribu si
et seulement si elle vérifie les trois propriétés suivantes :
1. Ω ∈ A.
11
12 CHAPITRE 1. RAPPELS DE MATHÉMATIQUES ET PROBABILITÉ
– A0 la tribu engendrée la famille des singletons {ω} de Ω. Cette tribu est utile lors de
la détermination d’une loi de probabilité.
– La tribu complète P(Ω).
Dans le cas où Ω est fini ou infini dénombrable (ce qui sera toujours le cas dans ce qui
suit), alors ces deux tribus sont identiques.
Définition 1.1.2 (Espace probabilisable). Le couple (Ω, A) est appelé espace probabi-
lisable. Dans le cas où Ω est fini dénombrable, A = P(Ω).
1.1.3 Probabilités...
Maintenant que l’on a défini la structure de l’univers dans lequel se déroule l’expérience
aléatoire qui nous intéresse, reste à donner une forme au hasard. C’est ce qu’on appelle
probabiliser l’espace probabilisable. Il s’agit simplement de définir la probabilité qu’un
événement A ∈ A survienne.
Définition 1.1.3 (Probabilité). P est une probabilité définie sur l’espace probabilisable
(Ω, A) si et seulement si P est une application de A vers R qui vérifie les propriétés
suivantes :
1. 0 ≤ P (A) ≤ 1, ∀A ∈ A.
2. P (Ω) = 1 (Axiome de normalisation).
3. Pour toute famille finie (Ai )0≤i≤n d’événements de A, deux à deux incompatibles,
on a :
n n
!
[ X
P Ai = P (Ai )
i=1 i=1
4. Pour toute famille dénombrable (Ai )i∈N d’événements de A, deux à deux incom-
patibles, on a :
∞ ∞
!
[ X
P Ai = P (Ai )
i=1 i=1
(Axiome de σ-additivité).
Le nombre réel du second membre est la somme de la série de terme général P (Ai ).
Dans la mesure où les termes de cette série sont positifs et que la probabilité de l’union
de n Ai est majorée par 1, cette série est toujours convergente.
Définition 1.1.5. On sait que P (Ω) = 1, mais on peut trouver des événements A 6= Ω
et tels que P (A) = 1. On dit que ces événements sont quasi-certains ou presque sûrs.
On sait que P ( ) = 0, mais on peut trouver des événements A 6= et tels que P (A) =
0. On dit que ces événements sont quasi-impossibles ou négligeables.
Définition 1.1.6. Deux distributions P et P 0 sont dites équivalentes si elles ont les
mêmes négligeables.
Définition 1.1.7 (Espace probabilisé). Le couple (Ω, A, P ) est appelé espace probabi-
lisé.
Notons finalement que la donnée d’une probabilité P sur un espace probabilisable est
équivalent à la donnée d’une distribution de probabilité sur Ω.
X:Ω→R
où fx est la densité de probabilité de la variable aléatoire X. Cette densité est telle
que :
Z i+
P (i− < X < i+ ) = f (x)dx (1.1)
i−
On reviendra plus tard sur la définition d’une densité. Le moment d’ordre 1 est l’espérance :
Z
E [X] = xfx (x)dx (1.2)
Ω
La covariance est une mesure de dépendance entre deux variables aléatoires. Soit deux
variables aléatoires X et Y :
Z Z
Cov(X, Y ) = (x − E[x]) (y − E[y]) f (x, y)dx, dy (1.6)
X(ω) Y (ω)
1.1. DES VARIABLES ALÉATOIRES ET DES HOMMES 15
où f (x, y) est la densité de la loi jointe de X et Y . Elle mesure la chance qu’on deux
séries d’évoluer de concert. Il est aisé d’interpréter son signe, mais pas son amplitude.
En normant la covariance par le produit des écart-types de X et Y , on obtient une
mesure dont il est aisé d’interpréter la valeur : le coefficient de corrélation. Il se calcule
comme suit :
Cov(X, Y )
ρ(X, Y ) = (1.7)
σX σY
ρ(X, Y ) ∈ [−1; 1], ce qui rend son interprétation aisée.
∀a ∈ R, F (a) = P (X ≤ a)
où f est une densité de X. En effet, tout autre fonction g de R dans R, qui coincide
avec f sauf sur un ensemble fini de points de R est aussi une densité de probabilité de X.
On ne propose pas revue des principales distributions, dans la mesure où il est aisé de
trouver ces informations sur Wikipédia.
16 CHAPITRE 1. RAPPELS DE MATHÉMATIQUES ET PROBABILITÉ
Supposons que l’on ait affaire à un groupe d’individu composé d’hommes et de femmes,
de bons et de mauvais élèves. On peut s’interesser à la probabilité de tirer au hasard un
bon éléve au sein de cette population, mais on peut aussi s’intéresser au fait de tirer un
bon éléve parmi les hommes. Il s’agit ici de la probabilité de tirer un bon éléve, sachant
que l’on tire parmi les hommes. On parle dans ce cas de probabilité conditionnelle. On
note :
La règle de Bayes permet de faire le lien entre les probabilités conditionnelles et non
conditionnelles :
T
P (A B)
Définition 1.1.12 (Probabilité conditionnelle). P (A|B) = P (B) .
P (B|A)∗P (A)
Définition 1.1.13 (Règle de Bayes). P (A|B) = P (B)
Dans le cas continu, il est possible d’obtenir une densité conditionnelle. Soit un couple
de variables aléatoires (X, Y ) définies sur un espace probabilisé (Ω, A, P ). Alors la
densité de X sachant Y s’écrit :
fX,Y
fX|Y =
fY
A ceci s’ajoute une propriété importante : la loi des espérances itérées.
Définition 1.1.14 (Espérances itérées). Soit une variable aléatoire X sur un espace
probabilisé. Alors sont espérance peut être calculée comme suit :
où Y est une autre variable aléatoire par rapport à laquelle on conditionne.
Théorème 1.1.1. Soit une variable aléatoire continue X, ayant fX (.) pour densité et
soit le support de X suivant :
X = {x|fX (x) ≥ 0} (1.14)
Si h(.) est une fonction dérivable et strictement monotone de domaine X et d’image
U, alors U = h(X) a pour densité :
dx
fU (u) = fX (h−1 (u)) , ∀u ∈ U (1.15)
du
= 0 sinon (1.16)
Il existe enfin une version réelle de cette fonction caractéristique que l’on appelle fonc-
tion génératrice des moments.
Définition 1.1.16. La fonction caractéristique d’une variable aléatoire X est la fonc-
tion suivante :
φ(t) = E[etX ] (1.20)
On l’appelle fonction génératrice des moments du fait de la propriété suivante :
Proposition 1.1.2. Soit une variable aléatoire X de fonction génératrice des moments
φ(t). Alors on a :
∂ k φ(t)
E[X k ] = (1.21)
∂tk t=0
Exercice : déterminer les deux premiers moments d’une loi normale à partir de la
fonction génératrice des moments..
18 CHAPITRE 1. RAPPELS DE MATHÉMATIQUES ET PROBABILITÉ
Théorème 1.2.3 (Théorème central limite). Soit (Xi )i∈N une suite de variables aléatoires
indépendantes et identiquement distribuées, avec ∀i, E[Xi ] = m, V[Xi ] = σ 2 . On a
alors :
1 Pn
n i=1 Xi − m L
σ −→ N (0, 1) (1.28)
√
n
Une matrice carrée est telle que n = p. Une matrice symétrique est une matrice carrée
telle que mi,j = mj,i , ∀i 6= j.
Le rang (colonne) d’une matrice M(n×p) est le nombre maximum de colonnes qui sont
linéairement indépendantes les unes des autres. En notant r(A) le rang d’une matrice
A qui soit une M(n × p), il vient naturellement que :
Une matrice carré d’ordre n est non singulière si son rang est égal à n (par exemple
une matrice diagonale).
Définition 1.3.1 (Matrice inverse). Soit A une matrice carrée d’ordre n. Son inverse,
notée A−1 si elle existe, est la matrice de même dimension telle que :
Si la matrice A−1 existe, alors on dit que la matrice A est inversible. Cette matrice
existe si et seulement si la matrice A est plein rang, autrement dit si la matrice A est
non singulière.
20 CHAPITRE 1. RAPPELS DE MATHÉMATIQUES ET PROBABILITÉ
Dans ce qui suit, on suppose acquit les éléments suivants : la somme de deux matrices,
le produit de deux matrices, la trace d’une matrice ainsi que le déterminant d’une ma-
trice. On rappelle en revanche différentes opérations de différenciation de matrices.
y = Ax (1.32)
y = xT A (1.34)
alors on a :
∂y
= AT . (1.35)
∂x
Soit maintenant :
y = xT Ax (1.36)
On dit que y est la variable expliquée ou endogène et {x1 , x2 , ..., xn } sont les variables
explicatives ou exogènes. est une perturbation aléatoire : il vient perturber une re-
lation qui, sans lui, resterait stable. Ce terme reçoit de nombreuses dénominations,
selon les champs d’application de l’économétrie ainsi que les séries étudiées. Quelques
exemples : il est possible de qualifier les de bruit (artefact statistique qui ne comporte
pas d’information particulière), d’erreur de mesure (erreur sur la compréhension de y
que permet de le modèle), de choc exogène (un choc qui ne transite pas par les variables
du modèle)...
Afin d’estimer cette relation, on utilise un échantillon de l’ensemble des variables. On
note yi la ième valeur de l’échantillon des y et xi,j la ième valeur de l’échantillon de la
jème variable. La valeur observée de yi est alors la somme de deux composantes : l’une
déterministe (les xi,j , ∀j) et l’autre aléatoire, i .
L’objectif est d’estimer les paramètres inconnus du modèle, d’utiliser les données pour
étudier la validité de certaines propositions théoriques, et éventuellement de former une
21
22CHAPITRE 2. RETOUR SUR LE MODÈLE LINÉAIRE : CAS UNIVARIÉ ET MULTIVARIÉ
prévision de y. On s’appuie dans cette démarche sur un corpus statistique bien connu :
les estimateurs du maximum de vraisemblance.
H 2 (Plein rang). Il n’existe pas de relation linéaire entre les variables indépendantes.
Si X est la matrice des observations, on aura notamment X 0 X de plein rang, et donc
inversible.
H 3 (Exogéneité des variables indépendantes). E[i |xi,1 , xi,2 , ..., xi,n ] = 0. L’espérance
de la perturbation conditionnellement aux variables exogènes est nulle : les variables
exogènes n’apportent plus aucune information sur la perturbation.
Une fois ces hypothèses mises à jour, on revient sur l’écriture du modèle. On preferera
travailler avec des matrices, plutôt qu’avec un indiçage couteux en place et en patience
(la mienne). Mieux, la plupart des logiciels de statistique/économétrie ont le bon goût
de fonctionner également en matriciel. Cette démarche simplifie grandement les calculs,
comme on le verra par la suite.
Soit x:k le vecteur colonne de T observations de la variable xk , k = 1, .., n. Soit X un
matrice MT,n constituée par la concaténation des différents vecteurs colonnes. Dans la
plupart des cas, la première colonne de X est constituée par un vecteur colonne unitaire
(1, 1, ..., 1)0 , de façon à ce que β1 soit la constante du modèle.
Sur la base de ces éléments, il est possible de réécrire 2.2 sous forme matricielle :
X β 0 + |{z}
y = |{z} (2.4)
|{z} |{z}
MT,1 MT,n Mn,1 MT,1
– L’hypothèse H2 rappelle qu’il ne peut exister de relation linéaire entre les variables
explicatives. X est une matrice MT,n : elle doit donc être de rang n, i.e. de plein
rang colonne. Deux conditions sont à remplir pour cela :
– une condition d’identification : il est nécessaire de disposer de n observations au
moins ;
– la non-colinéarité entre vecteurs colonne.
Rajoutons également qu’au moins un régresseur doit varier (et par conséquent être
non constant). Dans le cas contraire, la condition de plein rang n’est pas vérfiée (deux
colonnes sont identiques à une constante mutliplicative près).
L’hypothèse de nullité de l’espérance des erreurs n’est pas une hypothèse contraigante
(voir Greene (2002), page 15).
E[1 1 |X] E[1 2 |X] ... E[1 n |X]
E[2 1 |X] E[2 2 |X] ... E[2 n |X]
E[0 |X] = (2.9)
.. .. .. ..
. . . .
E[n 1 |X] E[n 2 |X] ... E[n n |X]
2
σ 0 ... 0
0 σ 2 ... 0
= . (2.10)
. .. .. ..
. . . .
0 0 ... σ 2
1 0 ... 0
0 1 ... 0
= σ2 . . . . (2.11)
.
. . . . . . .
0 0 ... 1
= σ2I (2.12)
– Notons enfin que les perturbations sont supposées suivre une loi normale d’espérance
nulle et de variance égale à σ 2 I. Il s’agit d’une hypothèse bien fondée étant donné la
structure de (les perturbations sont formées d’une suite de chocs, de même loi et
de mêmes moments 1 et 2 : le théorème central limit s’applique sans restriction).
M in (Y − X βb0 )2 (2.15)
La résolution est simple. Il suffit de dériver l’expression à minimiser par rapport à β et
de chercher la valeur de β (unique avec nos conditions) l’annulant :
2X 0 (Y − X βb0 ) = 0 (2.16)
⇔X 0 Y = X 0 X βb0 (2.17)
⇔βb0 = (X 0 X)−1 X 0 Y (2.18)
L’estimateur des paramètres du modèle par la méthode MCO est alors :
Théorème 2.1.1 (Régression orthogonale). Si les variables dans une régression mut-
liple ne sont pas corrélées (autrement dit, si elles sont orthogonales), alors les esti-
mations obtenues sont les mêmes que celles obtenues dans les regressions individuelles
simples.
Quelle est alors la solution pour β2 ? Les équations normales (i.e. l’équation obtenue
après dérivation de l’erreur quadratique telle qu’elle est définie par le modèle) sont alors
les suivantes :
Ce problème peut être résolu de deux façons possibles : soit en utilisant les règles
connues sur les matrices partionnées, soit en développant l’expression matricielle.
A−1 −1
−A−1
A11 A12 11 (I + A12 F2 A21 A11 11 A12 F2
= −1 (2.23)
A21 A22 −F2 A21 A11 F2
Avec :
X1 (X10 X1 )−1 X10 = P1 : il s’agit du projecteur dans l’espace des X1 . Il s’agit d’une
matrice symétrique et idempotente [i.e. P 2 = P ]. On remarque :
y = X1 β10 + (2.31)
= X1 (X10 X1 )−1 X10 y + (2.32)
= P1 y + (2.33)
⇔ y − = P1 y (2.34)
Théorème 2.1.4 (Gauss Markov). Sous les hypothèses du modèle linéaire, l’estima-
teur des moindres carrés ordinaire d’un modèle linéaire est optimal dans la classe des
estimateurs sans biais, conditionnellement aux régresseurs.
SCE
R2 = (2.36)
SCT
SCT − SCR
= (2.37)
SCT
Pn
− Ȳ )2 − ni=1 (Yi − Ybi )2
P
i=1 (Yi P
= n 2
(2.38)
i=1 (Yi − Ȳi )
Dans le précédent calcul, SCT représente la somme des carrés totaux, SCR, la somme
des carrés résiduels et SCE, la différence entre les deux, c’est à dire la somme des
carrés expliqués. Ȳ est la moyenne empirique de Y . Le R2 se définit donc comme le
rapport de la somme des carrés expliqués sur la somme des carrés totaux. L’idée est en
fait de décomposer la variance de Y en une variance expliquée par le modèle estimé et
une variance qui n’a pas pu être expliquée. Naturellement, plus R2 est grand et plus -
en principe - le modèle peut être soupçonné d’être explicatif de la variable endogène.
Cet indicateur, par construction est toujours compris entre 0 et 1. Ainsi plus le R2 est
proche de 1 (de 0) et plus (moins) le modèle est explicatif de Y .
Une mise en garde s’impose ici : dans une regression multiple, le R2 augmente natu-
rellement lorsqu’une variable supplémentaire (vérifiant certaines conditions qui ne sont
pas détaillées ici) est ajoutée. Ceci signifie donc que l’intoduction d’un grand nombre
de variable peut naturellement conduire à obtenir un R2 important, quand bien même
le pouvoir explicatif du modèle est médiocre.
R2 N − P − 1
Ftest = ∼ F (P, N − P − 1) (2.39)
1 − R2 P
On en déduit aisément :
βb
p ∼ Tn−p−1 (2.41)
(X 0 X)−1 σ
Pour n − p − 1 grand (supérieur à 30), on a Tn−p−1 → N (0, 1). Là encore, lorsque la
valeur de la statistique de test est supérieure à la valeur critique pour un niveau de
risque défini, on rejette l’hypothèse nulle de nullité du coefficient du paramètre. Ainsi,
si le quantile de la loi de Student à 95% et n − p − 1 degrés de liberté est plus petit que
la valeur de la statistique calculée ( √ 0β −1 ), on est conduit à rejeter l’hypothèse
b
(X X) σ
de nullité du paramètre. Le paramètre est alors significativement différent de 0. En
pratique, on compare la valeur de cette statistique à 2 : il s’agit de la valeur la plus
courante du quantile d’une loi de Student. La règle est donc : si la valeur de la statistique
calculée est inférieure à 2, alors αi = 0 ; dans le cas contraire, αi = αbi . L’intuition est
donc : on conduit un test de Student pour être bien sûr que la valeur estimée par MCO
soit bien différente de 0. Il s’agit de vérifier la qualité de l’estimation.
Nota Bene : 1. Intuitivement, plus la variance des résidus σ est importante (autrement
dit, moins le modèle semble être explicatif du comportement de Y ) et plus l’erreur
possible lors de l’estimation des paramètres est potentiellement importante.
Pn
i − b
i=1 (b i−1 )2
d= Pn 2 (2.42)
i=1 bi
d ∼ 2(1 − ρ) (2.43)
Dans les deux derniers cas, l’estimation par les moindres carrés ordinaires n’est pas sa-
tisfaisante. Il est alors nécessaire de développer des méthodes plus avancées permettant
d’intégrer l’existence de cette autocorrélation.
T
1X
µk = (Xi − X̄)k (2.44)
T
i=1
La skweness s’estime donc par µ3 /µ32 et la kurtosis par µ4 /µ22 . La statistique de Jarque
et Berra vaut donc :
T 2 T
s= Sk + (Ku − 3)2 → χ2 (2) (2.45)
6 24
Le qqplot
Il est également possible de procéder à un test d’adéquation des résidus à une loi
paramétrique quelconque. On se reportera au chapitre 1 pour la méthodologie.
30CHAPITRE 2. RETOUR SUR LE MODÈLE LINÉAIRE : CAS UNIVARIÉ ET MULTIVARIÉ
n
Y
f (y1 , ..., yn |θ) = f (yi |θ) = L(θ, y) (2.46)
i=1
Cette densité jointe est la fonction de vraisemblance, définie comme une fonction du
vecteur de paramètres inconnus (θ), où y indique les données observées (qui ne sont
donc pas une inconnue). On remarque que l’on note la densité jointe comme une fonc-
tion des données conditionnellement aux paramètres alors que, lorsque l’on forme la
fonction de vraisemblance, on note cette fonction en sens inverse, comme une fonction
de paramètres conditionnellement aux données observées. Dans ce qui suit, on suppose
que les paramètres sont constants et inconnus : l’enjeux de la méthode est d’exploi-
ter l’information disponible dans l’échantillon afin d’en inférer une valeur probable des
paramètres.
Il est généralement plus simple de travailler avec le logarithme de la fonction de vrai-
semblance :
n
X
lnL(θ|y) = lnf (yi |θ) (2.47)
i=1
On parle dans ce cas de log-vraisemblance. Ajoutons qu’il est courant de travailler sur
la densité d’un processus conditionnellement à un autre processus. C’est du moins ce
qui se passe dans le modèle linéraire : les erreurs sont bien i.i.d., ce qui fait que y|x est
aussi un processus iid. Soit le modèle linéaire gaussien suivant :
y = Xβ 0 + (2.48)
⇔yi = β1 + β2 x1,i + ... + βp xp−1,i + i (2.49)
n
X
lnL(θ|y,X) = lnf (yi |Xi , θ) (2.50)
i=1
n
(yi − x:,i β 0 )2
1X 2
=− lnσ + ln(2π) + (2.51)
2 σ2
i=1
n
1 X (yi − x:,i β 0 )2
1 2
= − nlnσ − nln(2π) − (2.52)
2 2 σ2
i=1
e−θ θyi
f (yi |θ) = (2.53)
yi !
Puisque les observations sont i.i.d., leur densité jointe, qui est la vraisemblance de cet
échantillon, est :
1 P1
Y e−10θ θ i=1 0yi e−10θ θ20
f (y1 , ..., y1 0|θ) = 0f (yi |θ) = Q1 = (2.54)
i=1 0yi !
207, 360
i=1
n
X n
X
lnL(θ|y) = −nθ + lnθ yi − ln(yi !) (2.55)
i=1 i=1
n
∂lnL(θ|y) 1 X
= −n + yi = 0 ⇔ θ̂EM V = ȳn (2.56)
∂θ θ
i=1
les mêmes que celles qui maximisent lnL(θ|y). La condition nécessaire pour maximiser
lnL(θ|y) est :
∂lnL(θ|y)
=0 (2.57)
∂θ
Il s’agit de l’équation de vraisemblance. Le résultat général est que l’EMV est une racine
de l’équation de vraisemblance. L’application aux paramètres du processus générant les
données d’une variable aléatoire discrète suggèrent que le maximum de vraisemblance
est une bonne utilisation des données. Cette intuition reste à généraliser dans ce qui
suit.
yi = xi β 0 + i (2.58)
n 2
n 1X xi − µ
lnL(θ|x) = −nln(σ) − ln(2π) − (2.59)
2 2 σ
i=1
Dans le cas de , on a :
n
n 1 X i 2
lnL(θ|x) = −nln(σ) − ln(2π) − (2.60)
2 2 σ
i=1
Il est alors nécessaire de savoir passer de la loi de à celle de y, lorsque l’on a spécifié
le modèle. Pour cela, on utilise un changement de variable qui est rappelé ici :
Proposition 2.2.2 (Changement de variable). Si x est une variable aléatoire de fonc-
tion de répartition fx (.) et si y = g(x), alors la densité de y s’exprime comme suit :
Cette proposition sera particulièrement importante pour l’analyse des séries temporelles
(troisième partie de cet opus). Ici, la transformation de à y est = y − Xβ 0 , donc la
∂
jacobienne (matrice des dérivées premières) est égale à l’unité ( ∂y = 1). On en déduit
naturellement la vraisemblance associée à y :
n 2
yi − xi β 0
n 1X
lnL(θ|x) = −nln(σ) − ln(2π) − (2.62)
2 2 σ
i=1
1 Pn
" #
∂lnL xi (y − xi β 0 )
∂β σ2 i=1 P 0
∂lnL = n
(y −x β 0 )2 = (2.63)
∂σ − 2σn2 + 12 i=1 σi4 i 0
Pn Pn
∂ 2 lnL ∂ 2 lnL x2i
" # " #
i=1 xi i
∂β∂β 0 ∂β∂σ −P i=1
σ2
− σ 4n
∂ 2 lnL ∂ 2 lnL = n P 2 (2.66)
i=1 xi i n i=1 i
∂σ∂β 0 ∂σ∂σ − σ4 2σ 4
− σ6
X0X
σ2
0
I(θ) = n (2.67)
0 2σ 4
σ 2 (X 0 X)−1
−1 0
V[θ] = I(θ) = 2σ 4 (2.68)
0 n
LR → χ2 (J) (2.70)
y 0 = x0 β 0 + 0 (2.71)
Par le théorème de Gauss Markov, on a :
e0 = y 0 − ŷ 0 (2.73)
= x0 (β 0 − β̂ 0 ) + 0 (2.74)
Où β est la vraie valeur du paramètre et β̂ son estimateur. On en déduit alors aisément
la variance de la prédiction :
q
0 0 0
ICα = ŷ ± tα/2 V[e |X, x ]
b (2.79)
Afin de déterminer la qualité d’une prévision, on utilise en général les deux statistiques
suivantes :
s
1 X
RM SE = (yi − ŷi )2 (2.80)
n0
i
1 X
M AE = 0 |yi − ŷi | (2.81)
n
i
2.4. UNE CALIBRATION SIMPLE DU CAPM 37
On parle de Root Mean Square Error et de Mean Absolute Error. n0 désigne le nombre
de périodes de prévisions.
On cherche donc à minimiser la somme des erreurs commises pour chacune des obser-
vations de la rentabilité du titre i, élevées au carré. On cherche donc :
X
M in 2i (2.86)
X
⇔ (ri − rf − β(rm − rf ))2 (2.87)
Pour trouver un minimum sur une fonction convexe, il suffit d’égaler la dérivée à 0. Ici,
on :
∂SSR(β) X
= (rm − rf )(ri − rf − βi (rm − rf )) (2.88)
∂β
X
= β(rm − rf )2 − (ri − rf )(rm − rf ) (2.89)
X X
= β(rm − rf )2 − (ri − rf )(rm − rf ) (2.90)
38CHAPITRE 2. RETOUR SUR LE MODÈLE LINÉAIRE : CAS UNIVARIÉ ET MULTIVARIÉ
En notant r̃i et r̃m les rentabilités du titre i ainsi que du marché dont on a retrancher
le taux sans risque, il vient une écriture simple de l’estimateur des MCO, qui peut se
réécrire de façon matricielle aisément :
P
r̃i r̃m
β= P (2.94)
(r̃m )2
En notant Ri la matrice n×1 contenant les n observations de rentabilités du titre (dont
on a retranché le taux sans risque) i et Rm de même taille, celle contenant celles du
marché, il est alors possible de réécrire ces sommes de façon matricielle (faire les calculs
pour s’en convaincre) :
T
β = (Rm Rm )−1 RmT Ri (2.95)
Cet estimateur est l’estimateur MCO de β. On admettra qu’il possède la distribution
suivante :
T
βb ∼ N (β, (Rm Rm )−1 σ2 ) (2.96)
σ2 est la variance du terme d’erreur et β la vraie valeur du paramètre à estimer. On
est ainsi en mesure de construire un test de Student permettant de tester sur l’estima-
tion de β est différente de 0 avec une probabilité importante. On se borne ici à fournir
la méthodologie générale permettant de construire ce test : il s’agit pour nous d’une
recette de cuisine financière dont les fondements ne sont pas démontrés.
β̂
Pour conduire ce test, il suffit de comparer à 1,96 (le quantile à 95% d’une loi
V[β̂]
normale). Si la valeur est supérieure, alors β̂ 6= 0.
2.4.5 Le R2
Le R2 ou coefficient de détermination est une mesure de la qualité globale du modèle
proposé. Il s’agit du rapport entre la variance de l’estimation de l’on donne du modèle
et la variance de la variable expliquée. Il s’agit donc du rapport :
V[βi rm ]
R2 = (2.107)
V[ri ]
Dans le cas du MEDAF, il est aisé de calculer ce R2 à partir du β :
βi2 σm2
R2 = = ρi,m (2.108)
σi2
On aboutit donc à une expression simple de ce R2 .
40CHAPITRE 2. RETOUR SUR LE MODÈLE LINÉAIRE : CAS UNIVARIÉ ET MULTIVARIÉ
# Dans ce programme, on estime les beta du CAPM par MCO. On teste l’existence de
# alpha. Enfin, on estime la SML par MCO.
capm<-function(x,Rf){
# Formattage de la base de données
x=[Link](x)
[Link]=x
x=x-Rf
SP=cbind(matrix(1,nrow(x),1),x[,1])
titres=x[,2:ncol(x)]
xx=solve(t(SP)%*%SP)
theta=xx%*%(t(SP)%*%titres)
res=titres-SP%*%theta
[Link]=apply(res,2,var)
test=cbind(diag(xx)[1]*[Link]([Link]),diag(xx)[2]*[Link]([Link]))
test=sqrt(test)
test=t(theta)/test
ptest=pnorm(test)
#Calcul de la SML
[Link]=t(theta)[,2]
[Link]=[Link](apply([Link][,2:ncol([Link])],2,mean))
pente=sum([Link]*[Link])/sum([Link]^2)
[Link]=([Link])
[Link]=[Link]*pente
par(bg="lightyellow")
plot([Link],[Link],type="l",col="red",ylab="Rentabilité",xlab="Beta",
main="Estimation de la SML")
lines([Link],[Link],type="p",col="blue")
# Calcul des R2
R2=[Link]([Link]^2)*var(SP[,2])/apply(titres,2,var)
return(list(theta=theta, R2=R2, test=test, ptest=ptest))
}
2.4. UNE CALIBRATION SIMPLE DU CAPM 41
Estimation de la SML
●
0.020
●
●
● ●
0.018
●
● ●
●
●
●
0.016
●
Rentabilité
● ●
● ●
0.014
●
●
● ●
● ●
0.012
●
0.010
●
●
●
●
Beta
On propose dans ce qui suit quelques extensions du modèle linéraire standard : les
modèles non linéaires et les modèles à équations multiples. L’exposé s’appuie à la fois sur
Greene (2002) et Davidson and MacKinnon (1993). Il est à noter qu’il existe une version
française de ce dernier ouvrage librement accessible sur le site de Russel Davidson
(voir l’url fournie en bibliographie). Il s’agit d’une introduction succinte : tout lecteur
soucieux de dépasser ce stade se reportera avec profit aux deux références qui sont
fournies plus haut.
y = x(β) + (3.1)
avec x(β) une forme fonctionnelle liant les variables exogènes x et le vecteur des pa-
ramètres β. est une perturbation i.i.d. de moyenne nulle et de variance égale à σ 2 .
La fonction scalaire x(β) est une fonction de regression (non linéaire) qui détermine
l’espérance de y conditionnellement à β et à x. Ici encore, un modèle de régression non
linéaire doit être identifié si l’on désire obtenir des estimations uniques des paramètres.
Classiquement, l’estimation peut se faire par moindres carrés. La fonction objectif est
alors :
n
X
M in (yi − xi (β))2 (3.2)
i=1
43
44 CHAPITRE 3. EXTENSIONS DU MODÈLE DE BASE
où y désigne une matrice M(n × 1) composé des différentes valeurs de yi et x(β) une
matrice M(n × p) composé des n fonctions de regression xi (β), i = 1, ...p. Cette somme
des carrés peut être explicité comme suit :
∂x(β) 0 ∂x(β) 0
−2y +2 x(β) = 0 (3.5)
∂β ∂β
∂x(β) 0
(y − x(β)) = 0 (3.6)
∂β
On retrouve ici la condition d’orthogonalité entre les résidus et la dérivée des régresseurs,
une version modifiée de la condition d’orthogonalité entre regresseurs et résidus dans le
cas MCO. Le problème est qu’il n’existe pas formule analytique permettant d’obtenir
l’expression des estimateurs comme dans le cas MCO : il est alors nécessaire d’utiliser
des algorithmes d’optimisation permettant de trouver le minimum du programme men-
tionné plus haut. Nous reviendrons sur ces algorithme plus loin : il seront également
utiles pour l’estimation des modèles de séries temporelles.
Un dernier point reste à noter dans ce bref exposé des moindres carrés non linéaires :
les conditions de premier ordre sont nécessaire, mais pas suffisante pour garantir le fait
que β̂ soit un minimum. Il peut exister plusieurs valeurs de β qui vérifient les conditions
de premier ordre, mais qui soient des minima lcaux, des points stationnaires ou même
des maxima locaux. En définitive, rien ne garantit que la fonction à minimiser soit
globalement convexe, i.e. :
∂ 2 SSR(β)
Hij (β) = (3.7)
∂βi βj
Rappel 2. Une matrice M M(n×n) est dite définie positive si ∀x, une matrice colonne
composée de n éléments réels on a :
x0 M x > 0 (3.8)
Une autre façon de prouver le caractère défini positif d’une matrice est de montrer
que l’ensemble de ses valeurs propres est strictement positif. Enfin, si l’ensemble des
sous matrices d’une matrice carré admet un déterminant positif, alors cette matrice est
définie positive. Il s’agit de la généralisation matricielle de la postivité d’un scalaire.
Une matrice définie négative est une matrice dont l’opposé est définie positive.
3.1. MODÈLE DE RÉGRESSION NON LINÉAIRE 45
Cette condition n’est assurée que pour quelques cas particuliers, dont les MCO :
∂ 2 SSR(β)
Hij (β) = = −x0 x (3.9)
∂βi βj
Il est aisé de prouver d’après ce qui vient d’être rappelé que cette matrice est bien
définie positive.
Terminon enfin par une courte discussion des conditions d’identification de ce type
de modèle. On distingue deux srotes d’identification, l’identification locale et l’iden-
tification globale. Les estimation des moindres carrés non linéaire ne seront identifiés
localement qu’à la condition que pour tout modification infinitésimale de β̂, la valeur
de la fonction objectif s’élève. Ceci revient à supposer que la matrice hessienne est
strictement convexe en β, de sorte que :
Pour une petite variation δ. Cette condition est analogue au caractère défini positif
de la matrice hessienne. La stricte convexité implique que SSR(β) soit incurvée dans
toutes les directions ; aucun plat n’est autorisé quelle que soit la direction. Si SSR(β)
était plate dans une direction donnée au voisinage de β̂, il serait possible de s’éloigner
de β̂ dans cette direction, sans jamais modifier la valeur de la somme des résidus au
carré (du fait des conditions du premier ordre). Par conséquent, β̂ ne sera pas l’unique
estimateur des moindres carrés non linéaires.
L’identification locale n’est cependant pas suffisante. Une condition plus générale est
l’identification globale :
y = βγ + γ 2 xt + (3.12)
y = AK β1 Lβ2 (3.13)
46 CHAPITRE 3. EXTENSIONS DU MODÈLE DE BASE
qui redevient un modèle linéaire dans le cas où l’on passe au log :
y1 = X1 β10 + 1 (3.15)
y2 = X2 β20 + 2 (3.16)
... (3.17)
0
ym = Xm βm + m (3.18)
On suppose que n observations sont utilisées pour l’estimation des paramètres des m
équations. Chaque équation a Km régresseurs, pour un total de K = m
P
K
i=1 j . On pose
T > Ki . On suppose que les perturbations ne sont pas corrélées entre observations. En
conséquence,
σ11 In σ12 In ... σ1m In
σ21 In σ22 In ... σ2m In
E[0 |X1 , ..., Xm ] = Ω = = Σ ⊗ In (3.25)
..
.
σm1 In σm2 In ... σmm In
Chaque équation est une régression classique. Les paramètres peuvent donc être estimés
de manière convergent par la méthode MCO. La régression généralisée s’applique aux
données dites empilées :
y11
.. β1,1
.
β1,2
y1n X1 0n,k2 ... 0n,km
..
.
y21 0n,k X2 ... 0n,km
1
.. = .. β1,k1 (3.26)
.. ..
. . . ... .
β2,1
y2n 0n,k 0 n,k ... Xm
..
1 2
.
..
.
βm,km
ymn
1,1
1,2
..
.
1,n
+ (3.27)
2,1
..
.
m,n
Y = Xβ 0 + (3.28)
Où Y est une matrice M(n×m, 1), X une matrice M(n×m, K), β 0 une matrice M(K, 1)
et une matrice M(n × m, 1). Une fois ce travail préliminaire accompli, tournons nous
vers l’estimation de ce type de modèle.
A ce stade plusieurs stratégies d’inférence sont envisageables. On montre dans ce qui suit
que l’estimation par maximum de vraisemblance ou par Moindres Carrés Généralisés
(MCG), comme dans le cas simple des MCO, coincident exactement. On s’intéresse
finalement au cas où Σ est bloc diagonal.
(Y 0 − βX 0 )Σ−1 ⊗ In (Y − Xβ 0 ) (3.30)
0 −1 0 −1 0 −1 0 0 −1 0
=Y Σ ⊗ In Y − βX Σ ⊗ In Y − Y Σ ⊗ In Xβ + βX Σ ⊗ In Xβ (3.31)
En dérivant la précédente expression, on obtient :
Il est aisé de s’en convaincre sur un exemple de faible dimensions, par exemple pour
m = 2, K = 2. Autre remarque : en spécifiant Σ = σIm , on retrouve exactement
l’estimateur des moindres carrés. On détaille quelques propriétés des estimateurs MCG :
Proposition 3.2.2 (Estimateur sans biais). L’estimateur MCG est un estimateur sans
biais de β.
Pour finir, on a supposé que jusque ici que Σ était connue. Ceci n’est généralement pas
le cas. On procède donc à une estimation préliminaire de Σ : les résidus des moindres
carrés peuvent être utilisées pour estimer la matrice de variance-covariance des résidus.
En note ˆ les résidus des MCO, l’estimation de Σ est alors :
1 0
Σ̂ = ˆ ˆ (3.37)
n
On parle alors d’estimateur des moindres carrés quasi-généralisés pour β̂ estimé par
MCG, une fois l’estimation de Σ accomplie. Les propriétés asymptotiques de l’estima-
teur MCQG ne proviennent pas de l’estimateur sans biais de Σ ; seule la convergence
est nécessaire. L’estimateur a les mêmes propriétés que l’estimateur MCG.
– Si les équations ont les mêmes variables explicatives, alors MCO et MCG sont iden-
tiques.
Sur la base des notations introduites plus haut, la vraisemblance concentrée du modèle
SURE s’écrit :
n 1
lnL = − ln(|Ω|) − (Y − Xβ 0 )0 Ω−1 (Y − Xβ 0 ) (3.39)
2 2
En dérivant par rapport à β 0 , il vient :
X 0 Ω−1 (Y − Xβ 0 ) = 0 (3.40)
0 −1 0 −1 0
⇔X Ω Y =XΩ Xβ (3.41)
0 0 −1 −1 0 −1
⇔β = (X Ω X) XΩ Y (3.42)
On retrouve bien l’estimateur des MCG : celui profite donc de l’ensemble des propriétés
des estimateurs du maximum, developpées plus haut. L’estimateur de la variance se
déduit des propriétés suivantes :
Proposition 3.2.5 (Dérivation matricielle 1). Soit A une matrice carrée, inversible et
de déterminant positif. Alors :
∂ln(|A|)
= (A0 )−1 (3.43)
A
Proposition 3.2.6 (Dérivation matricielle 2).
Une fois que l’on connait ces propriétés, il est aisé de voir que :
∂lnL n 1
−1
= Ω − (Y − Xβ 0 )(Y − Xβ 0 )0 (3.45)
∂Ω 2 2
En égalisant la dernière dérivée à 0, il vient :
1 0
Ω= (3.46)
n
où = Y − Xβ 0 . On retrouve bien un estimateur habituel de la variance.
On commence donc par réécrire l’ensemble des matrices de façon appropriée. Avec
l’estimation de Σ, on fournit une estimation de Ω :
Ω = Σ ⊗ In (3.47)
Il ne reste alors plus qu’à estimer la matrice β, qui comporte (si tout se passe bien)
n × m éléments en colonne. L’estimateur βb est alors :
La fonction capm.R a été complétée pour prendre en compte l’estimation par MCQG : le
code est fournit ci-après. Une mise en garde s’impose : le calcul de Ω conduit à construire
une matrice aux dimensions imposantes. Il est par conséquent possible que de nombreux
PC ne puisse pas permettre l’estimation par MCQG d’un système d’équations où m est
grand. Dans notre cas, pour 30 titres et 250 dates, la matrice Ω est de taille 7500×7500 :
elle comporte donc... 56 250 000 éléments ! Il est par conséquent nécessaire de disposer
d’une mémoire vive... très importante.
Les resultats retournés par la procédure sont fournis dans la table 3.2.4 et comparés à
ceux obtenus par MCO. Le constat principal est la suivant : les MCQG n’apporte rien
en comparaison des MCO lors de l’estimation du CAPM. Ceci tient à la matrice de
variance-covariance des erreurs MCO qui est une matrice diagonale presque surement.
Pour tester la valeur d’une corrélation ρ, on rappelle que la statistique de test est :
√ ρ
Tρ = n − 2p (3.49)
1 − ρ2
Cette statistique suit, sous H0 : ρ = 0, une loi de student à n − 2 degrés de liberté (où
n est le nombre de données disponibles). On compare donc Tρ à 1,96, comme on l’a fait
pour les tests de Student.
Finalement, ceci ne sert qu’à montrer que dans de nombreux cas pour lesquels la struc-
ture de corrélation se réduit aux simples variances, il n’est pas nécessairement utile de
recourir aux MCG/MCQG : les MCO suffisent amplement. Pour refaire ces estimations
vous-mêmes, il suffit d’ajouter à la fonction capm.R les quelques lignes qui suivent, ainsi
que de modifier return, comme c’est ici le cas :
for (i in 1:ncol(titres)){
[Link][(n*(i-1)+1):(i*n),(2*(i-1)+1):(2*i)]=SP
}
# Calcul de la matrice de variance covariance des résidus MCO
[Link]=var(res)
omega=solve(kronecker([Link],diag(1,nrow(titres),nrow(titres))))
# Estimation des paramètres
[Link]=t(solve(t([Link])%*%omega%*%[Link])%*%(t([Link])%*%omega%*%[Link]))
return(list(theta=theta, R2=R2, test=test, ptest=ptest, [Link]=[Link]))
Alcoa AT.T Boeing Caterpilar Chevron [Link] Disney DuPont [Link]
MCO alpha 0.003024439 -0.01330993 0.004927275 -0.001995679 -0.006321155 -0.001599203 0.0084873 0.001150471 -0.01899075
beta 1.029791342 0.74692300 1.143201992 1.054370309 0.840468354 0.802655308 1.1226112 1.016821780 0.59013381
MCQG alpha 0.003024439 -0.01330993 0.004927275 -0.001995679 -0.006321155 -0.001599203 0.0084873 0.001150471 -0.01899075
beta 1.029791 0.746923 1.143202 1.054370 0.8404684 0.8026553 1.122611 1.016822 0.5901338
Optimisation de fonctions à
plusieurs variables par algorithme
Ne le cachons pas, proposer ce type de chapitre n’est pas monnaie courante : rares sont
les cours de statistique à insister sur de tels aspects computationnels. Je suis intime-
ment convaincu que de trop nombreux économètres ne programment pas leurs propres
procédures d’estimation, et ceci est extrêment dommageable. Des générations d’élèves
sont formés dans l’idée que l’économétrie se fait simplement en lançant des procédures
SAS ou pire : Eviews. Il existe un certain illétrisme économétrique parmis beaucoup
d’élèves : il est possible de programmer ses propres procédures d’estimation (notamment
en R) et ce sans faire d’efforts incroyables. Je propose dans ce qui suit les principales
intuitions nécessaires à la programmation de fonction permettant d’optimiser un fonc-
tion de p variables. Bien évidemment, on ne propose ici pas l’ombre d’une preuve des
méthodes proposées : il s’agit bien plutôt d’une chapitre de cuisine économétrique. Il
existe de nombreuses références proposant les preuves de ce qui va suivre, mais je doute
qu’un élève de M1 en tire un quelconque profit.
Les différentes recettes proposées ont été tirées pelle-melle de : Harvey (1990), Des-
champs (2004) Quinn (2001) et Greene (2002). Notez que ce chapitre du Green est li-
brement téléchargeable sur le site de Pearson Education (et en Français !)1 . On propose
trois types de méthodes : une première approche intuitive vise à déterminer graphique-
ment (méthode de recherche par quadrillage) les fonctions de un ou deux variables. Ce
type de méthode laisse ensuite la place à l’ensemble des méthodes basées sur le gra-
dient : plus grande pente, Newton Raphson, méthode du score et BHHH. Ces méthodes
se heurtent cependant aux problèmes liées aux multimodalités du maximum de vrai-
semblance : les méthodes de type recuit simulé sont brièvement présentées.
1
[Link]
55
56CHAPITRE 4. OPTIMISATION DE FONCTIONS À PLUSIEURS VARIABLES PAR ALGORITHME
C’est ce qui est fait par le code suivant (à ceci pret que la fonction proposée représente
également la logvraisemblance en 3D).
grid<-function(x,n){
mu<-seq(1,3,length=n)
sigma<-seq(1,3,length=n)
T=nrow(x);
lnL=matrix(0,n,n);
for (i in 1:n){
for (j in 1:n){
mm2=sum((x-mu[i])^2);
4.1. POUR COMMENCER... 57
lnL[i,j]=-T*log(sigma[j])-1/2*(mm2/(sigma[j]^2));
}}
persp(mu,sigma,lnL, theta = 50, phi = 30, expand = 0.5, col = "lightblue");
max=min(lnL);
for (i in 1:n){
for (j in 1:n){
if(max<lnL[i,j]){max=lnL[i,j]; indexi=i;indexj=j}
}}
mumax=mu[indexi];
sigmamax=sigma[indexj];
return(list(mu=mu,sigma=sigma,lnL=lnL,sol=cbind(max,mumax,sigmamax)))
}
On obtient alors les graphiques présentés en figure 4.1 et 4.2, en utilisant les fonction
persp et contour de R. On en déduit alors l’estimation de nos paramètres. L’idée est
essentiellement de se construire une représentation mentale de ce qu’est une vraisem-
blance (ça existe !), avant de passer à des méthodes plus avancées.
Logvraisemblance en 3D
lnL
a
u
m
sig
3.0
2.5
2.0
1.5
1.0
θt+1 = θt + λt ∆t (4.3)
4.2. LES MÉTHODES DU GRADIENT 59
On développe ici quelques méthodes bien connues : méthode de la plus grande pente,
méthode de Newton Raphson et méthode du score par BHHH. L’ensemble de ces
méthodes sont dires méthodes du gradient dans la mesure où ∆t est systématiquement
de la forme :
∆ t = W t Gt (4.4)
−G0 G
λt = (4.6)
G0t Ht Gt
où H est la matrice hessienne, i.e. la matrice des dérivées secondes de f . On a donc :
∂2f
H= (4.7)
∂θ∂θ0
Ainsi l’itération de la plus grande pente est alors :
G0 G
θt+1 = θt − Gt (4.8)
G0t Ht Gt
En supposant que l’on a p paramètres à estimer, θ est alors une M(p, 1), G est également
une M(p, 1) et H est une M(p, p). On vérifie donc que les dimensions des matrices cor-
respondent bien.
1. Etablir un θ0
2. Déterminer G0 et H0
3. Déterminer ∆t et λt
60CHAPITRE 4. OPTIMISATION DE FONCTIONS À PLUSIEURS VARIABLES PAR ALGORITHME
4. Déterminer θ1
P 2
5. Tester G < . Si oui, stop. Si non, on continue la boucle...
steep<-function(theta,x,iter){
H=matrix(0,2,2)
G=matrix(1,2,1)
n=nrow(x)
check=matrix(0,2,1)
i=1;
while(sum(G^2)>0.0000001){
mm=sum(x-theta[1,1])
mm2=sum((x-theta[1,1])^2)
H[1,1]=-n/(theta[2,1]^2);
H[2,2]=H[1,1]-3*(mm2)/(theta[2,1]^4);
H[1,2]=-(2*mm)/(theta[2,1]^3);
H[2,1]=H[1,2];
G[1,1]=mm/(theta[2,1]^2);
G[2,1]=-n/theta[2,1]+(mm2)/(theta[2,1]^3);
cat(i,"\n");
check=cbind(check,theta-solve(H)%*%G);
theta=[Link]((t(G)%*%G)/(t(G)%*%H%*%G))*G;
i=i+1
}
return(list(theta=theta,check=check))
}
Ici, on a fixé à 0.0000001, de façon arbitraire. On illustre les résultats de cette méthode
sur notre problème initial de vraisemblance gaussienne en figure 4.3. On a fixé volon-
tairement θ0 loin de la vraie valeur des paramètres : on observe une certain nombre de
sauts, même si on trouve au final la bonne valeur approchée des paramètres.
1000
0
d$ch[1, ]
−1000
−2000
0 10 20 30 40
Index
0 = G + Hh (4.11)
−1
⇔h = −H G (4.12)
En d’autres termes, la direction optimale pour optimiser f (retour sur notre problème)
est de choisir θt+1 tel que :
1. Etablir un θ0
2. Déterminer G0 et H0
3. Déterminer ∆t et λt
4. Déterminer θ1
P 2
5. Tester G < . Si oui, stop. Si non, on continue la boucle...
newton<-function(theta,x,iter){
H=matrix(0,2,2)
G=matrix(1,2,1)
n=nrow(x)
62CHAPITRE 4. OPTIMISATION DE FONCTIONS À PLUSIEURS VARIABLES PAR ALGORITHME
check=theta
i=1;
while(sum(G^2)>0.0001){
mm=sum(x-theta[1,1])
mm2=sum((x-theta[1,1])^2)
H[1,1]=-n/(theta[2,1]^2);
H[2,2]=H[1,1]-3*(mm2)/(theta[2,1]^4);
H[1,2]=-(2*mm)/(theta[2,1]^3);
H[2,1]=H[1,2];
G[1,1]=mm/(theta[2,1]^2);
G[2,1]=-n/theta[2,1]+(mm2)/(theta[2,1]^3);
cat(i,"\n");
check=cbind(check,theta-solve(H)%*%G);
theta=theta-solve(H)%*%G;
i=i+1
}
return(list(theta=theta,check=check))
}
On applique également cette méthode à notre problème, en partant encore d’un point
éloigné du véritable maximum. Là, l’algorithme de Newton Raphson diverge et aboutit
à la solution m = −1182571474 et σ = −168165462, soit des valeurs abhérantes. Ceci
est représenté sur la figure 4.4.
−8.0e+08
−1.2e+09
5 10 15 20
Index
Cette méthode converge cependant dans de nombreux cas. Si la fonction est quadra-
tique, elle atteint l’optimum en une itération depuis n’importe quel point de départ.
Si la fonction est globalement concave, cette méthode reste probablement la meilleure.
Elle est particulièrement adaptée à l’estimation du maximum de vraisemblance.
4.2. LES MÉTHODES DU GRADIENT 63
1. Etablir un θ0
2. Déterminer G0 et BHHH0
3. Déterminer ∆t et λt
4. Déterminer θ1
P 2
5. Tester G < . Si oui, stop. Si non, on continue la boucle...
HHH<-function(theta,x){
G=matrix(1,2,1)
n=nrow(x)
check=theta
i=1;
while(sum(G^2)>0.0001){
mm=sum(x-theta[1,1])
mm2=sum((x-theta[1,1])^2)
BHHH=cbind((x-theta[1,1])/(theta[2,1]^2),-1/theta[2,1]+((x-theta[1,1])^2)/(theta[2,1]^3
H=(t(BHHH)%*%BHHH);
G[1,1]=mm/(theta[2,1]^2);
G[2,1]=-n/theta[2,1]+(mm2)/(theta[2,1]^3);
cat(theta,"\n");
check=cbind(check,theta+solve(H)%*%G);
theta=theta+solve(H)%*%G;
i=i+1
}
return(list(theta=theta,check=check))
}
64CHAPITRE 4. OPTIMISATION DE FONCTIONS À PLUSIEURS VARIABLES PAR ALGORITHME
5
0
2 4 6 8 10 12 14
Index
1. Choisir θ0
2. Tirer autant de que de paramètres
3. Déterminer θ1
4. Comparer lnL(θ0 ) et lnL(θ1 )
5. Conserver le θi rendant la log-vraisemblance maximale
6. Recommencer la boucle
random<-function(theta,x,iter){
nb=length(theta)
logv<-function(x,theta){
mm2=sum((x-theta[1,1])^2);
lnL=-nrow(x)*log(theta[2,1])-1/2*(mm2/theta[2,1]^2)
return(list(lnL=lnL))
}
thetachain=theta
thetanew=matrix(0,2,1)
for(i in 1:iter){
for (j in 1:nb){thetanew[j,1]=theta[j,1]+rnorm(1)/2}
if (logv(x,thetanew)$lnL>logv(x,theta)$lnL){theta=thetanew}
cat(i,"\n")
thetachain=cbind(thetachain,theta)
}
return(list(theta=theta,lnL=logv(x,theta)$lnL,thetachain=thetachain))
}
Appliqué à notre problème, ce code fournit une estimation correcte de nos paramètres.
En partant d’un point éloigné de maximum, l’algorithme converge cependant lentement.
C’est ce qu’on observe sur la figure 4.6.
Il est possible de converger plus rapidement, en effectuant des saut dans les paramètres
plus importants. Cependant, ce faisant, on diminue également la probabilité de trouver
un point qui maximise la vraisemblance. Notons que l’algorithme s’alourdit considérablement
en présence d’un grand nombre de paramètres. En effet, on travaille à chaque itération
autour du voisinage du dernier maximum trouvé par l’algorithme : un voisinage à n
paramètres est bien plus grand qu’un voisinage à deux paramètres. Pour finir, il est
possible de coupler les méthodes basées sur le gradient (définissant la direction) avec un
algorithme construit pour générer des nombres positifs aléatoirement : on génère alors
une taille de pas aléatoire. L’idéal est de générer au début de l’algorithme des taille de
pas importantes, et de diminuer ces tailles au fur et à mesure de l’algorithme. Ceci est
d’ailleurs proche de ce qu’on appelle la température dans l’Algorithme de Métropolis
Hastings.
66CHAPITRE 4. OPTIMISATION DE FONCTIONS À PLUSIEURS VARIABLES PAR ALGORITHME
100
80
d$thetach[1, ]
60
40
20
0
Index
∗
θt+1 = θt + t (4.18)
∗
Ce qui change par rapport à la précédente méthode est le mode de selection entre θt+1
et θt . On compare la log-vraisemblance associée à chacun d’eux :
∗ ) > lnL(θ ), alors on retient θ
– si lnL(θt+1 ∗
t t+1 = θt+1 .
∗
– sinon, on ne rejette pas nécessairement θt+1 . On l’accepte avec une probabilité de
0lnL(θ∗ )−lnL(θ )
la forme : min(exp(− 1
T ), 1). Dans la pratique, on procède à un tirage
d’une loi uniforme, puis on compare ce tirage avec le précédent calcul. Si le tirage est
∗ , sinon, on conserve θ . La température T du système permet
inférieur, on choisit θt+1 t
d’augmenter cette probabilité d’accepter θt+1∗ . En général, il s’agit d’une fonction du
1. Déterminer θ0
2. Tirer et déterminer θ1∗
3. Comparaison de la lnL associée aux deux paramètres
4. Si lnL(θ1∗ ) > lnL(θ0 ) ⇒ θ1 = θ1∗
5. Sinon : on tire u selon une loi uniforme.
lnL(θ1∗ )−lnL(θ0 )
6. On calcule min(exp(− T ), 1) et on compare à u
lnL(θ1∗ )−lnL(θ0 )
7. Si u < min(exp(− T ), 1), alors θ1 = θ1∗ .
8. Sinon θ1 = θ0
9. On recommence la boucle...
randommh<-function(theta,x,iter){
nb=length(theta)
logv<-function(x,theta){
mm2=sum((x-theta[1,1])^2);
lnL=-nrow(x)*log(theta[2,1])-1/2*(mm2/theta[2,1]^2)
return(list(lnL=lnL))
}
thetachain=theta
thetanew=matrix(0,2,1)
for(i in 1:iter){thetanew[2,1]=-1;
while(thetanew[2,1]<0.1){for (j in 1:nb){thetanew[j,1]=theta[j,1]+rnorm(1)}}
if (logv(x,thetanew)$lnL>logv(x,theta)$lnL){theta=thetanew; cat("-","\n")}
if (logv(x,thetanew)$lnL<logv(x,theta)$lnL){u=runif(1);
if(u<min(exp((logv(x,thetanew)$lnL-logv(x,theta)$lnL))*exp(1+i^(1/100)),1)){theta=theta
}
cat(i,"\n")
thetachain=cbind(thetachain,theta);
plot(thetachain[1,],type="l",col="blue")
}
return(list(theta=theta,lnL=logv(x,theta)$lnL,thetachain=thetachain))
}
Il existe des façons bien plus complexe de contrôler la température du système que celle
proposée ici qui est déterministe. Graphiquement, sur un problème aussi simple que le
notre, il n’y a pas grande différence entre la précédente méthode et celle-ci (cf. figure
4.7).
68CHAPITRE 4. OPTIMISATION DE FONCTIONS À PLUSIEURS VARIABLES PAR ALGORITHME
30
20
10
Index
xt est appelé variable aléatoire. En principe, il n’existe pas ou peu de différence entre
les séries temporelles et l’économétrie, sinon que les variables sont indicées par t plutôt
que par i. Par exemple, si yt est généré par :
alors une estimation par MCO permet d’obtenir des estimateurs consistants, de la même
façon que dans le cas classique (i.e. indexé par i).
Heureusement pour nous, il est très rare que des séries temporelles soient i.i.d. (indépendantes
et identiquement distribuées) : il s’agit précisement de ce qui les rend intéressantes. Le
PIB en donnée trimestrielle constitue un bon exemple de cette dépendance temporelle :
en général, quand le PIB est anormalement haut en t (i.e. au-dessus de sa moyenne
historique ou non-conditionnelle), il y a de très fortes chances pour que la prochaine
valeur de ce même PIB soit elle aussi anormalement haute.
L’idée est donc de parvenir à caractériser la loi jointe de {..., xt−1 , xt , xt+1 , ...}. Evi-
dement, il serait interessant de mettre en oeuvre des méthodes non paramétriques
69
70CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Dans ce qui suit, on présente deux classes de modèles : l’une s’attachant à modéliser la
moyenne conditionnelle des processus (les modèles ARM A) et l’autre s’intéressant à la
modélisation de la variance conditionnelle (les modèles ARCH − GARCH).
Notons néanmoins que la théorie financière repose néanmoins sur une hypothèse proche
du bruit blanc. Le modèle de Black-Scholes repose sur l’hypothèse d’une diffusion sui-
vant un brownien géométrique pour le prix des actifs. On rappelle que si St est le cours
du sous-jacent dans le modèle de Black-Scholes, alors sa diffusion est de la forme :
On en déduit aisément :
Z t+1 t+1 t+1
σ2
Z Z
dlog(Ss ) = (µ − )ds + σdWs (5.7)
t t 2 t
2
St+1 σ
⇔log =µ− + σ(Wt+1 − Wt ) (5.8)
St 2
Sachant que (Wt+1 − Wt ) ∼ N (0, 1) (l’intervalle de temps est l’unité), on retrouve donc
bien un processus pour les rendements qui est à peu de choses pret un bruit blanc.
La seule différence entre les deux processus tient au drift obtenu après application de
la formule d’Îto. Malheureusement, le cours des actifs est rarement bruit blanc, ce qui
ne va pas sans poser quelques problèmes lors de l’utilisation de la formule de Black
et Scholes. Elle constitue néanmoins un benchmark intéressant pour l’évaluation des
options.
Comme on peut le constater, il s’agit à chaque fois d’une recette à base de bruits blancs
passés et de valeurs de départ pour xt . L’ensemble de ces modèles ont une moyenne
nulle, et sont utilisés pour représenter l’écart des séries à leur moyenne. Par exemple,
si une série {xt } a une moyenne égale à x̄ et suit un processus AR(1), alors :
est équivalent à :
Ainsi, la constante c absorbe l’effet moyen. On travaille dans ce qui suit principalement
à l’aide de modèle excluant la constante : elle est aisément estimable.
72CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
De façon plus formelle, l’opérateur retard est un opérateur qui produit une nouvelle série
de donnée (retardée) à partir d’une série {xt }. A partir de cette définition sommaire,
il est aisé de voir que :
On a donc :
En utilisant ces notations, il est alors possible de réécrire des modèles ARMA comme
suit :
(5.25)
ou plus simplement :
AR : a(L)xt = t (5.26)
MA : xt = b(L)t (5.27)
ARMA(p,q) : a(L)xt = b(L)t (5.28)
– une représentation à l’aide d’un polynome retard le plus petit possible est toujours
plus aisé ;
– les modèles AR sont les plus aisés à estimer par MCO ;
– les MA représentent xt en fonction de variables indépendantes : dans de nombreux
cas, ceci facilitera les calculs de variance et de covariance, comme on le verra plus loin.
5.2. LES MODÈLES ARMA 73
On retrouve donc le résultat précédent. On a supposé que |φL| < 1. Tous les processus
ARMA n’ont pas de représentation inversible (on parle d’inversibilité des processus ou
de processus inversible) de xt en fonction du passé de t .
Il est possible d’étendre les résultats, en montrant sous quelle condition un AR(p) peut
admettre une représentation MA(∞). Les calculs sont cependant un peu plus longs et
fastidieux. Un lecteur intéressé lira avec intérêt Cochrane (2005), page 14 et suivantes.
Dans un cas où E[xt ] = 0, ∀t, alors γj = cov(xt , xt−j ) = E[xt xt−j ]. On a de plus
γ0 = V[xt ].
On en déduit aisément l’expression du coefficient de corrélation :
γj γj
ρj = = (5.38)
γ0 V[xt ]
γ0 = σ2 (5.39)
γj = 0, ∀j 6= 0 (5.40)
ρ0 = 1 (5.41)
ρj = 0, ∀j 6= 0 (5.42)
On observe ainsi que les autocorrélations d’un MA(1) s’annulent à partir de l’ordre 1.
C’est ce qu’on observe sur les figure 5.2. Notons que la figure 5.1 présente l’allure d’un
processus bruit blanc, MA(1), AR(1) et ARMA(1,1).
5.2. LES MODÈLES ARMA 75
ρ0 = 1 (5.54)
θ1 + θ1 θ 2
ρ1 = (5.55)
1 + θ12 + θ22
θ2
ρ2 = (5.56)
1 + θ12 + θ22
ρi = 0, ∀i > 2 (5.57)
On a alors :
q q
" # !
X X
i
γ0 = V[xt ] = V (θi L )t = θj2 σ2 (5.59)
i=0 i=0
q q q
" #
X X X
i i
γk = E[xt xt−k ] = E (θi L )t (θi L )t−k = θi θi+k σ2 , ∀k ≤ q (5.60)
i=0 i=0 i=0
γk = 0, ∀k > q (5.61)
Remarque 2. Il y a une leçon importante à retenir de tout ceci : les calculs de co-
variance pour les processus MA est simple dans la mesure où les termes en E[j k ]
deviennent rapidement nuls quand j et k sont éloignés. Il n’en va pas de même des
processus AR.
φk
γk = σ 2 , ρk = φk (5.65)
1 − φ2
L’autre façon de retrouver ces résultats est de travailler directement sur xt , sans utiliser
l’astuce de l’inversion. Pour le même modèle AR(1) que précédement, on a :
Ainsi l’ACF d’un modèle AR est général un mélange entre une sinusoide et une expo-
nentielle, selon les signes de coefficient du modèle AR. Si les signes sont négatifs, les
autocorrélations paires seront positives, et celles impaires seront négatives. Quoi qu’il
arrive, |φ| < 1, d’où limk→0 φk = 0. Une décroissance relativement lente vers zéro se
traduit par un convergence en forme d’exponentiele de l’ACF vers 0.
Ceci tient donc pour m = 0. Dans un tel cas, p(k) ∈ [−1; 1]. On ajoute la propriété
suivante :
|Pk∗ |
p(k) = (5.73)
|Pk |
avec
1 ρ1 . . . ρk−1
ρ1 1 . . . ρk−2
Pk = (5.74)
.. .. .. ..
. . . .
ρk−1 ... ... 1
et
1 ρ1 . . . ρ1
ρ1 1 . . . ρ2
Pk = (5.75)
.. .. .. ..
. . . .
ρk−1 . . . . . . ρk
L’idée est donc de caractériser la dépendance entre les différents retards en éliminant
l’impact à chaque fois des retards intermédiaires. La solution la plus simple consiste a
utilisée la définition avancée plus haut utilisant les moindres carrés de façon itérative,
pour les différents ordres de l’autocorrélation partielle. On en déduit (intuitivement)
que la PACF d’un AR(p) devrait être égale à 0 pour les autocorrélations partielles d’un
ordre supérieur à p.
Sans en faire la preuve formelle, on retrouve (théoriquement) dans le cas d’un processus
MA(q) un fonction d’autocorrélation partielle qui est un mélange d’une fonction expo-
nentielle et d’une sinusoı̈de. Il s’agit donc d’un cas symétrique de l’ACF pour les AR(p).
Cov[x
d t , xt−k ]
ρ̂k = (5.78)
Cov[x
d t , xt ]
D’après le théorème central limite, la variable centrée tρk suit une loi normale centrée
réduite :
ρ̂k − ρk L
tρk = p → N (O, 1) (5.79)
V[ρ̂k ]
autocorrel<-function(x,lag,series){
# This function computes the ACF for any series of data.
#Dimension setting
n=nrow(x);
N=n-lag;
# Scaling of the data
x=scale(x)
#Creation of the matrix containing the lagged series
data=x[(lag+1):n,1];
for (i in 1:lag){
data=cbind(data,x[(lag-i+1):(n-i),1])
}
#Computation of the correlations from lag 1 to lag "lag"
correl=t(data[,1])%*%data[,1]
for (i in 1:lag){
correl=cbind(correl,t(data[,1])%*%data[,(i+1)])
}
5.2. LES MODÈLES ARMA 79
correl=t([Link](correl))
#Normalization to obtain the ACF
correl=correl/correl[1,1]
#Computation of the variance of the estimates
variance=matrix(1/N,nrow(correl),1)
for (i in 1:lag){variance[(i+1),1]=1/N*(1+2*sum(correl[1:i,1]^2))}
#Computation of the intervall around the 0 null hypothesis
intervalleup=1.96*sqrt(variance);
intervalledown=-1.96*sqrt(variance);
#Plotting the results : ACF and tests
#Definition of the plot windows
mindata=min(intervalleup,intervalledown,correl);
maxdata=max(intervalleup,intervalledown,correl);
#Plot
par(bg="lightyellow")
plot(correl,type="h",col="blue",ylim=c(mindata,maxdata),main="Autocorrelation
Function",xlab="Lag",ylab=series,bg="red")
lines(matrix(0,nrow(correl),1),col="black")
lines(intervalleup,type="l",col="red");
lines(intervalledown,type="l",col="red");
return(list(correl=correl, interval=cbind(intervalleup,intervalledown)))
}
pautocorrel<-function(x,lag,series){
# This function computes the PACF for any series of data.
#Dimension setting
n=nrow(x);
N=n-lag;
# Scaling of the data
x=scale(x)
#Creation of the matrix containing the lagged series
data=x[(lag+1):n,1];
for (i in 1:lag){
data=cbind(data,x[(lag-i+1):(n-i),1])
}
80CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
L’idée derrière tout ceci est de travailler sur des processus dont la moyenne et la variance
sont constante au cours du temps. Notez que si la loi du processus est gaussienne, il y
a équivalence entre les deux définitions de la stationnarité. Détaillons les conditions :
la première de ces conditions suppose l’existence ou la convergence du moment d’ordre
2 (la variance pour un processus centré) ; les deux conditions suivantes supposent que
l’espérance et la covariance du processus sont constantes au cours du temps. Il n’y a
donc ni rupture dans la moyenne, ni rupture dans la structure de dépendance au cours
du temps.
Dernier point dans cet errata théorique, le théorème de Wold est le théorème fonda-
mental de l’analyse des séries temporelles stationnaires.
Théorème 5.2.1 (Théorème de Wold). Tout processus stationnaire d’ordre deux (xt )
peut être représenté sous la forme :
∞
X
xt = ψj t−j + κt (5.87)
j=0
P∞
où les paramètres ψj satisfont ψ0 = 1, ψj ∈ R∀j ∈ N∗ et 2
j=0 ψj < ∞ et où t est un
bruit blanc i.i.d..
On dit que la somme des chocs passés correspond à la composante linéaire stochastique
de xt . Le terme κt désigne la composante linéaire déterministe telle que Cov[κt , t ] =
0, ∀j ∈ Z. Il s’agit par conséquent d’une formalisation de ce qui a été déjà été dit sur
la transformation d’un processus AR en un MA(∞). On n’en fournit pas la preuve.
Bruit blanc AR(1)
2
4
1
2
0
ar
wn
0
−1
−2
−2
−4
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Index Time
MA(1) ARMA(1,1)
3
4
2
1
2
0
0
ma
arma
−1
−2
−3
−4
0 50 100 150 200 250 300 0 50 100 150 200 250 300
82CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Time Time
0.8
0.8
ACF
ACF
0.4
0.4
5.2. LES MODÈLES ARMA
0.0
0.0
0 5 10 15 20 0 5 10 15 20
Lag Lag
MA(1) ARMA(1,1)
1.0
0.8
0.6
ACF
ACF
0.4
0.2
0.0
−0.2
0 5 10 15 20 0 5 10 15 20
83
Lag Lag
0.8
0.10
0.6
0.4
0.00
Partial ACF
Partial ACF
0.2
0.0
−0.10
5 10 15 20 5 10 15 20
Lag Lag
MA(1) ARMA(1,1)
0.8
0.1
0.6
0.0
0.4
−0.1
0.2
Partial ACF
Partial ACF
−0.3
−0.2
5 10 15 20 5 10 15 20
84CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Lag Lag
xt = (1 − φL)−1 (µ + t ) (5.90)
2
= µ + t + φ(µ + t−1 ) + φ (µ + t−2 ) + ... (5.91)
∞
X ∞
X
=µ φi + φi t−i (5.92)
i=0 i=0
∞
µ X
= + φi t−i (5.93)
1−φ
i=0
E[xt ] = 0 (5.95)
ce qui n’est bien sur pas surprenant. L’espérance conditionnelle à xt−1 est différente de
cette dernière expression :
Ainsi conditionellement à xt−1 , l’espérance d’un AR(1) est différente de l’espérance non
conditionnelle. Il est possible de dérouler les mêmes calculs pour la variance condition-
nelle et non conditionnelle :
∞
µ X
V[xt ] = V[ + φi t−i ] (5.97)
1−φ
i=0
σ2
= (5.98)
1 − φ2
Là encore, les moments conditionnels et non conditionnels ne coincident pas. L’estima-
tion de ce type de processus nécessite de travailler non plus sur la vraisemblance mais
sur la vraisemblance non-conditionnelle. Harvey (1990) en rappelle le principe général :
dans le cas de la vraisemblance non conditionnelle avec observations i.i.d., il est pos-
sible d’écrire la loi jointe du processus comme le produit des lois pour chacune des
observations, du fait de la propriété d’indépendance des observations. Ici, on travaille
en relachant cette hypothèse. On utilise alors le fait que conditionnellement à l’obser-
vation du passé, les observations sont i.i.d.. Pour se faire, on applique la règle de Bayes
rappelée dans le chapitre 1, de façon à reproduire la décomposition proposée dans le
cas de la vraisemblance. Dans le cas où l’on dispose de trois observations {x1 , x2 , x3 } il
est alors possible d’écrire :
lnL = ln(f (x1 )) + ln(f (x2 |x1 )) + ln(f (x3 |x2 , x1 )) (5.102)
Il est alors possible d’estimer les paramètres en utilisant les méthodes proposées au
chapitre 4. On présente la méthode dans le cadre d’un modèle AR(1). Si t ∼ N (0, σ 2 ),
alors, en utilisant les calculs précédents des moments conditionnels, le processus xt =
xt−1 + t suit une loi normale d’espérance conditionnelle φxt−1 et de variance σ2 . Sa
log-vraisemblance s’écrit alors comme suit :
n
n−1 1 X
lnL(x, φ, σ2 ) = − ln(2π) − (n − 1)ln(σ) − 2 (xt − φxt−1 )2 + ln(f (x1 ))
2 2σ
t=2
(5.103)
5.2. LES MODÈLES ARMA 87
Dans la plupart des cas, on néglige le terme ln(f (x1 )) : dans le cas d’échantillons
importants, son influence est minime. Ceci peut se réécrire sous forme concentrée et
matricielle :
1
lnL(x, φ, σ2 ) = −(n − 1)ln(σ ) − 2 (Xt − φXt−1 )0 (Xt − φXt−1 ) (5.104)
2σ
où Xt et Xt−1 sont des matrices M(n − 1 × 1) contenant les observations du processus
(xt )t∈Z . Les équations normales sont alors :
∂lnL 1 0
= 2 Xt−1 (Xt − φXt−1 ) = 0 (5.105)
∂φ σ
∂lnL n−1 1
=− + 3 (Xt − φXt−1 )0 (Xt − φXt−1 ) = 0 (5.106)
∂σ σ σ
(5.107)
Ce système admet la solution suivante :
0
θ = (Xt−1 Xt−1 )−1 Xt−1
0
Xt (5.108)
r
1 0
σ = t (5.109)
n−1 t
On retrouve donc exactement la même solution que dans le cas des MCO. Tout ce qui
a été dit précédement sur ces estimateurs est donc valable : on n’y reviendra pas.
où Φ est la matrice M(1 × p) des p coefficients à estimer et Xt−1:t−p est la matrice
M(n − p × p) des séries de données retardées.
Il est alors aisé de déterminer la log-vraisemblance conditionnelle et de retrouver les
estimateurs du maximum de vraisemblance, qui, là encore, coincident avec ceux des
MCO :
1
lnL(x, Φ, σ2 ) = −(n − 1)ln(σ ) − 2 (Xt − Xt−1:t−p Φ0 )0 (Xt − Xt−1:t−p Φ0 ) (5.114)
2σ
88CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
On commence comme précédement par constater que les moments d’ordre 1 et 2 condi-
tionnels et non conditionnels ne sont pas les mêmes :
Là encore les moments conditionnels et non conditionnels ne coincident pas. On rai-
sonnera là encore en terme de vraisemblance conditionnelle. On sait que :
[Link]<-function(n,theta){
epsilon=[Link](rnorm(n));
x=matrix(0,n,1);
for (i in 2:n){
x[i,1]=theta*epsilon[(i-1),1]+epsilon[i,1]
}
x[1,1]=epsilon[1,1]
return(list(x=x))
}
[Link]<-function(theta,x){
G=matrix(1,2,1)
n=nrow(x)
check=theta
i=1;
epsilon=matrix(0,n,1);
while(sum(G^2)>0.0001){
for (i in 2:n){epsilon[i,1]=x[i,1]-theta[1,1]*epsilon[i-1,1]}
BHHH=cbind((1/theta[2,1])*epsilon[1:(n-1),1]*(x[2:n,1]-theta[1,1]*epsilon[1:(n-1),1])
,-1/theta[2,1]+(1/theta[2,1]^3)*(x[2:n,1]-theta[1,1]*epsilon[1:(n-1),1])^2);
H=(t(BHHH)%*%BHHH);
G[1,1]=sum(BHHH[,1]);
G[2,1]=sum(BHHH[,2]);
cat(theta,"\n");
check=cbind(check,theta+solve(H)%*%G);
theta=theta+solve(H)%*%G;
i=i+1
}
plot(check[1,],type="l",col="blue",ylim=c(min(check),max(check)))
lines(check[2,],col="red")
return(list(theta=theta,check=check))
}
Enfin, la figure [Link] fournit la trajectoire pour différents points de départ des pa-
ramètres estimés.
1.0
0.05
0.8
0.00
0.6
MA(1) simulé
MA(1) simulé
−0.05
0.4
−0.10
0.2
0.0
−0.15
−0.2
5 10 15 20 5 10 15 20
Lag Lag
3
2
1
5.2. LES MODÈLES ARMA
x
0
−1
−2
−3
0 200 400 600 800 1000
Index
1.0
1.5
0.6
1.0
check[1, ]
check[1, ]
0.5
0.2
0.0
−0.2
5 10 15 0 10 20 30 40 50 60
Index Index
1.0
1.0
0.6
0.5
check[1, ]
check[1, ]
0.2
0.0
−0.2
−0.5
0 50 100 150 1 2 3 4 5 6 7
92CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Index Index
Là encore les t ne sont pas observables : il est nécessaire de procéder à leur estimation
à Θ = {θ1 , ..., θp } fixé. Il est alors possible d’écrire la vraisemblance associée à Θ ainsi
qu’aux observations. On détermine les t comme suit :
1. On suppose que les i pour i = {1, ..., q} sont nuls (espérance du processus bruit
blanc).
2. On détermine alors q+1 = xq+1 .
3. On peut alors déterminer q+2 = xq+2 − θ1 q+1 .
4. Puis q+3 = xq+3 − θ1 q+2 + θ2 q+1 .
5. On poursuit ainsi jusqu’à 2q+1 = x2q+1 − qi=1 θi 2q − i.
P
Une fois ceci fait, il est aisé de calculer la vraisemblance puis de la maximiser en utilisant
les méthodes proposées au chapitre 4. La loi conditionnelle de xt est :
Xq
xt |t−1 , ..., t−q ∼ N( θi t−i , σ2 ) (5.128)
i=1
Là encore, pour maximiser la vraisemblance, il suffit de mettre en oeuvre l’une des
méthodes présentées au cours du chapitre 4. On préfèrera naturellement utiliser des
méthodes de type scoring qui ont le bon gout de converger à tous les coups (en théorie).
L’utilisation de la matrice BHHH simplifie considérablement le calcul de la matrice
d’information de Fisher.
94CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Xp q
X
xt |xt−1 , ..., xt−p , t−1 , ..., t−q ∼ N( xt−i + t−i , σ2 ) (5.138)
i=1 i=1
Là encore, la difficulté tient au calcul des t−i liés à l’introduction d’une composante
MA(q). L’algorithme peut être de la forme :
Il suffit alors d’appliquer les méthodes présentées au chapitre 4 pour obtenir une esti-
mation des paramètres. Le code à developper est bien entendu plus complexe que ce
qui a été développé jusqu’ici. Les concepteurs de R proposent par défaut une fonction
permettant d’estimer les paramètres d’un modèle ARMA. Il s’agit de la fonction arima
qui utilise la synthaxe suivante1 :
Lorsque le processus est bien estimé, les résidus entre valeurs estimées et réelles par
le modèle doivent se comporter comme un bruit blanc. L’une des hypothèses de bruit
blanc est que l’espérance des résidus (et la moyenne empirique par conséquent) est
nulle. Si le processus (t ) est i.i.d., on doit alors :
n
1X
¯t = ˆt → 0 (5.143)
n
t=1
Les processus ARMA peuvent être vu comme une façon de stationnariser les séries,
abus de langage pour désigner le fait que l’on se débarrasse de la corrélation existant
dans les séries. Il est donc nécessaire de tester l’existence d’autocorrélation dans les
résidus. Si ce test conduit au constat que les résidus sont corrélés, c’est certainement
que le modèle est mal spécifié. Il existe différents tests d’autocorrélation :
– Test de Box et Pierce : l’idée de ce test est simple. Il est basé sur une somme des
carrés de autocorrélations pour un horizon allant de 1 à K. La statistique est la
suivante :
K
X
QBP = n ρ2k → χ2K−p−q (5.146)
k=1
5.2. LES MODÈLES ARMA 97
On rappelle ici le test de Jarque et Berra explicité plus haut. Celui-ci repose sur deux
statistiques connues pour la loi normale. La première est la skewness ou coefficient
d’asymétrie. Son expression est :
E[(X − E[X])3 ]
Sk = 3 (5.148)
σX
Cette statistique est une mesure de l’asymétrie de la distribution, i.e. de la façon dont
la densité s’étale de part et d’autre de son espérance. Dans le cas d’une loi normale,
cette statistique vaut 0 : la loi est parfaitement centrée.
La seconde statistique sur laquelle s’appuie le test de Jarque et Berra est la kurtosis
ou coefficient d’applatissement. Cette statistique mesure l’applatissement des queues
de distribution : plus celles-ci sont épaisses et moins le processus a de chances d’être
gaussien. Un loi normale a théoriquement une kurtosis égale à 3. On mesure cet index
à l’aide de la statistique suivante :
E[(X − E[X])4 ]
Ku = 4 (5.149)
σX
La méthode d’estimation utilisant une loi des erreurs possiblement différente de la vraie
loi de celle-ci, mais appartenant à la famille des lois exponentielles est appelée Pseudo
Maximum de Vraisemblance. L’estimation de la matrice de variance covariance des esti-
mateurs peut naturellement se faire en utilisant la matrice BHHH présentée au chapitre
précédent. Il s’agit d’une méthode particulièrement utile pour l’estimation des proces-
sus ARCH/GARCH que l’on verra par la suite.
1. La première démarche est celle exposée plus haut et fut la démarche fondatrice
de ces modèles, telles qu’elle fut proposée par Box et Jenkins. Il s’agit simplement
de sélectionner les ordres pour le processus AR et le processus MA au vu des au-
tocorrélogrammes du processus utilisé pour l’estimation. Il s’agit simplement du
prolongement naturel de ce qui a été dit lors de l’introduction du présent chapitre.
L’étude de l’ACF permet de sélectionner l’ordre du processus MA et l’étude la
PACF permet de sélectionner l’ordre du processus AR. Cependant, comme il l’a
été mentionné plus haut, un processus AR(1) présente une ACF avec une forte
persistance, ce qui n’est dû qu’à la ”contagion” de l’autoregressivité sur les er-
reurs du processus. Autrement dit, un AR(1) présente certes une ACF qui décroit
lentement, mais il s’agit simplement d’une persistance née de l’unique retard de
l’AR(1), perceptible sur la PACF. Moralité : si cette méthode fournie un point
de départ, elle est loin d’être suffisante.
φ̂i − φ
∼ N (0, 1) (5.151)
σφ
5.2. LES MODÈLES ARMA 99
Ce type de test a déjà fait l’objet de développements lors de l’exposé des méthodes
du maximum de vraisemblance.
– Le critère d’Akaike (AIC) : le meilleur des modèles ARMA est celui qui minimise
la statistique :
AIC(p, q) = T log(σ2 ) + 2(p + q) (5.152)
– Le critère bayésien (BIC) : le meilleur des modèles ARMA est celui qui minimise
la statistique :
2 p+q
BIC(p, q) =T log(σ ) − (n − p − q)log 1 − (5.153)
T
2
−1 σx
+ (p + q)log(T ) + log (p + q) (5.154)
σ2 − 1
– Le critère de Hanan et Quinn : le meilleur des modèles ARMA est celui qui
minimise la statistique :
2 log(T )
HQ(p, q) = T log(σ ) + (p + q)log (5.155)
T
Au final, lors de l’estimation de processus ARMA, l’ensemble de ces critères sont à
utiliser pour l’estimation des ordres p et q : on commence en général par l’étude de
l’ACF et de la PACF pour se faire un ordre d’idée sur le processus latent. On définit
en général un ordre maximum pour p et q, puis on procède de façon descendante : on
élimine progressivement (on parle d’approche stepwise) les paramètres non significatifs,
tout en conservant un oeil sur les critères d’information.
La simplicité de ces processus les a rendu très attrayants : il est possible de prévoir
n’importe quelle série autoregressive pour un ordre important, à partir de la seule
connaissance de son passé. Dans le cas d’un AR, l’estimation se fait très simplement,
par MCO. On montre encore qu’il est aisé d’obtenir un intervalle de confiance pour la
prévision, en utilisant le théorème de Wold.
Première remarque : non, le théorème de Wold ne sert pas à rien. Il permet de déterminer
l’intervalle de confiance de la prévision simplement en s’appuyant sur la représentation
MA des processus stationnaires. Deuxième remarque : l’erreur de prévision va gran-
dissante au fur et à mesure que l’on s’éloigne de la date t. Troisème remarque : il
est toujours possible d’obtenir cet intervalle de confiance en estimant un modèle MA
avec un ordre important afin d’obtenir les φi nécessaire à l’estimation de l’intervalle de
confiance.
D’après ce qui vient d’être dit, il est évident qu’il est possible de construire un intervalle
de confiance de la forme :
k−1
" #
X
ICα = x̂t+k ± tα/2 φ2i σ2 (5.174)
i=0
Malheureusement, qui dit simplicité dit également pauvreté de la prévision. Les modèles
ARMA sont incapables de percevoir un retournement de tendance, puisqu’on ne fait
que multiplier la tendance actuelle. Pour remédier à ceci, des modèles à seuil ont été
introduits. Il en sera question plus tard lors d’une simple application. Pire, les mouve-
ments de prix, de rentabilité ou bien encore les évolutions des chiffres de l’économie sont
le résultats d’une mutlitude de chocs, de dépendances conditionnnelles entre séries...
Bref, xt ne contient certainement pas assez d’information pour parvenir à construire
une prévision robuste de son avenir ! Là encore, pour remédier à ceci, des modèles mul-
tivariés ont été introduits. On traitera de ces modèles dans la section consacrée à la
prise en compte des liens entre variables macroéconomiques.
102CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
L’inflation est l’une des variables macroéconomiques les plus suivies par les marchés,
et notamment les marchés de taux. La raison à cela est simple : il s’agit d’une variable
suivie par les Banques Centrales afin d’établir leur politique. On parle d’instrument de
la politique monétaire : les Banques Centrales sont censées réagir par une montée des
taux à tout accroissement de l’inflation, de façon mécanique. En réalité, la véritable
variable cible des banquiers centraux est l’inflation qu’il anticipent à moyen terme.
Celle-ci n’est malheureusement pas mesurable : on a recourt en France et en Europe à
L’Indice des Prix à la Consommation ou IPC pour mesurer cette inflation.
La mesure d’inflation suivie par les Banques Centrales est en réalité le glissement men-
suel de l’indice des prix à la consommation. En notant Pt l’indice des prix à la consom-
mation en date t, le glissement mensuel se note comme suit :
Pt − Pt−12
πt = (5.176)
Pt−12
2.5
2.0
1.5
1.0
Time
Cette intuition est confirmée par l’étude des ACF/PACF qui exhibent naturellement
une persistance importante, pour des ordres de retard eux-mêmes importants. Ces gra-
phiques sont représentés en figure 5.8.
1.0
0.8
0.8
0.6
0.6
Partial ACF
0.4
ACF
0.4
0.2
0.2
0.0
0.0
−0.2
−0.2
Lag Lag
Conclusion de cette courte étude : l’indice des prix à la consommation semble se com-
porter comme une marche aléatoire. La meilleure prévision de l’inflation de demain
devrait ainsi être celle d’aujourd’hui. En oubliant la composante ma1, et en posant
φ = 1, on a alors :
πt = πt+1 + t (5.177)
où It est l’information disponible à la date t. Il s’agit d’un fait qui commence à être
bien connu dans le monde de l’économie monétaire.
3.0
2.5
2.0
IPC
1.5
1.0
Time
Quoiqu’il en soit, il n’est pas certain que l’inflation soit tout à fait une marche aléatoire
pour ce qui est de la zone euro depuis 2000. En effet, on remarque visuellement que
la moyenne pour la période 1997-2000 ne semble pas être la même que celle pour la
période 2000-2006 (cf figure 5.7). Confirmons cette intuition : étudions la structure de
la série de l’IPC depuis 2000. La série est présentée en figure 5.10.
2.0
Time
1.0
0.6
0.8
0.4
0.6
Partial ACF
0.2
0.4
ACF
0.2
0.0
0.0
−0.2
−0.2
0.0 0.5 1.0 1.5 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Lag Lag
Finalement, les résultats obtenus sont réellements différents des précédents : l’inflation
moyenne est plus importante que dans le cas précédent et le modèle est globalement un
modèle autoregressif qui n’est plus marche aléatoire (le terme en AR(2) vient ”corriger”
l’importance du 0,88. Au final, notre petite estimation souligne que la zone euro est ac-
tuellement dans un régime inflationniste plus important que son régime de long terme.
Ceci est certainement dû à l’effet Balassa-Samuelson : lors de la formation d’unions
monétaires, un effet rattrapage par le haut du niveau des prix se produit. L’intégration
de nouveaux pays à bas niveau de vie conduit à chaque fois à l’accroissement du niveau
général des prix en Europe (il s’agit d’un index prix en glissement annuel).
La figure 5.13 présente l’évolution (sur une base mensuelle) de ces taux depuis 2000.
Premier constat : il existe un nombre important de dates pour lesquelles le taux ne
change pas. En notant rt le taux cible à la date t, il existe de nombreuses dates pour
lesquelles on a :
rt = rt−1 (5.180)
106CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
2.0
Time
3.5
3.0
2.5
2.0
0 20 40 60 80
temps
Ceci risque bel-et-bien de faire apparaitre une correlation persistante dans l’ACF du
processus. On observe ceci en figure 5.14. En effet, on observe tout d’abord un phénomène
de contagion évident entre les différentes dates sur l’ACF qui décroit lentement vers
0. Au contraire, la PACF admet un unique pic, proche de 1, suggérant que l’on peut
écrire pour le taux BCE un modèle de la forme :
rt = φrt−1 + t (5.181)
– Il est tout d’abord possible d’estimer un modèle MA avec un lag important. Ceci
signifierai que chaque taux BCE peut s’écrire comme la somme des chocs passés sur
ces taux (une sorte de représentation MA(∞) d’un processus AR(1). Après estima-
tion, un modèle MA(15) semble convenir : il minimise le critère BIC. Les coefficients
estimés sont les suivantes :
Estimations T-stats
ma1 0,9115 6,80731889
ma2 0,8645 3,92954545
ma3 1,0792 5,96902655
ma4 1,4888 7,61145194
ma5 1,5917 6,49408405
ma6 1,783 4,97766611
ma7 1,5166 5,9732178
ma8 1,5969 6,67600334
ma9 1,6871 5,59568823
ma10 1,3082 4,19698428
ma11 0,9441 4,30506156
ma12 0,6917 3,56546392
ma13 1,0329 4,61116071
ma14 0,7159 3,71124935
ma15 0,2105 1,72258592
Constante 2,9544 12,8788143
– La seconde stratégie, qui est la plus naturelle, consiste à estimer un AR(1), tout
en se doutant que les taux suivent probablement un processus proche de la marche
aléatoire. L’estimation d’un simple AR(1) donne des résultats similaires au MA(15).
Les estimations sont :
AR(1) Constante
Estimations 0,9772 2,8363
T-stats 59,5853659 5,11229272
1.0
1.0
0.8
0.8
0.6
0.6
Partial ACF
0.4
ACF
0.4
0.2
0.2
0.0
0.0
−0.2
−0.2
0 5 10 15 5 10 15
Lag Lag
On présente en figure 5.15 les prix black scholes calibrés sur volatilité historique, im-
plicite globale et implicite locale.
Prix BS
100
80
60
option
40
20
0
0 50 100 150
Index
0.004
0.002
0.000
0 50 100 150
Index
0.6
0.8
0.6
0.4
Partial ACF
ACF
0.4
0.2
0.2
0.0
0.0
0 5 10 15 20 5 10 15 20
Lag Lag
Poussons tout de même l’analyse un peu plus loin. Les prix d’options génèrent en général
ce que l’on appelle un smile : la volatilité est une fonction non linéaire du strike. On
représente cette surface de volatilité en figure 5.18.
Volatilité implicite
sg[100
:120, 1
:8]
Tim
e to
ma
tur
ity
ik e
Str
Une piste d’explication possible pour rendre compte de cette relation est d’étudier les
composantes AR et MA de chacune des volatilité implicites. C’est ce qu’on représente
en figures 5.19 et figures 5.20 : on retrouve une relation non linéaire entre strike et
composante.
5.2. LES MODÈLES ARMA 111
0.90
0.88
0.86
−0.5
−0.6
−0.7
Cette liste est bien évidement non exhaustive. Elle a cependant le mérite de mettre en
lumière le fait que dès lors que l’on souhaite s’extraire des standards de la finance (la
marche aléatoire des théories de l’efficience et de l’évaluation d’option), il est nécessaire
de considérer l’ensemble de ces propriétés avec beaucoup d’attention.
Le fait que ce type d’index connaissent un succès évident ne fait que souligner le besoin
qu’on les marchés financiers de disposer d’intruments de mesure de la volatilité. Ce
index est l’une des tentatives d’établissement de ces mesures et semblent relativement
robuste (cf. revue de littérature de Poon (2005)). Il ne s’agit pas de la seule façon de
faire : on propose ici deux autres tentatives plus ou moins intéressantes : l’index high-
low et le carré des rendements.
(lnHt − lnLt )2
σt2 = (5.182)
4ln2
114CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Cet estimateur est celui proposé par Bollen and Inder (2002) dans le cadre de rendement
suivant un brownien géométrique. Il s’agit d’une application d’une mesure initiallement
proposée par Parkinson (1980). Garman and Klass (1980) propose une amélioration de
la mesure de Parkinson qui prend la forme suivante :
Ht 2 pt 2
2
σt = 0.5 ln − 0.39 ln (5.183)
Lt pt−1
Ces estimateurs de la volatilité sont assez sensibles aux valeurs extremes, c’est à dire
aux rendements anormalement importants (qu’il s’agissent de rendement négatifs ou
positifs). L’importance de la probabilité d’occurence des événements extrèmes rend le
processus effectivement suivi par les rendements incompatibles avec la loi normale. Dans
cette mesure, les mesures H-L constitue une approximation intéressante de la variance
des rendements.
Notons finalement que ces indices H-L modélisent d’une certaine façon la variance
conditionnelle du processus des rendements : celle-ci change chaque jour. Il n’est par
conséquent pas question d’utiliser ces résultats pour modéliser la variance de l’ensemble
de la série des rendements d’un titre donnée, ce qui n’aurait de toute façon pas grand
sens.
0.015
0.000
Time
Time
Index VIX
VIX measure
0.03
0.01
Time
0.4
0.0
Lag
0.4
0.0
Lag
Autocorrélations du VIX
0.8
ACF
0.4
0.0
Lag
Fig. 5.22 – ACF des différentes mesures de la volatilité appliquées au S&P 500
de la volatilité des rendements sur la base de données journalières. Par exemple, Lopez
(2001) propose le modèle suivant pour les rendements :
rt = µ + σt t (5.184)
Dans ce cas, le carré des rendements (en négligeant le drift) permet de mesurer la
volatilité de ces mêmes rendements. Cependant, 2t suit naturellement une loi du χ2 à
σ2
un degré de liberté, dont la médiane est 0,455 : il s’en suit que 2t est inférieur à 2t
dans plus de la moitié des cas. En effet, on a :
σt2 σt2
2 2 2
P rt < = P σt t < (5.186)
2 2
2 1
= P t < (5.187)
2
= 0.52 (5.188)
Autrement dit, le carré des rendements (en l’absence de drift) est un estimateur sans
biais de la variance des rendements. Cependant, dans plus de la moitié des cas, il sous-
estime la variance des résidus. Ce résultat semble s’affaiblir lorsqu’on utilise le carré
des rendements toutes les 5 minutes plutot que le carré des rendements journaliers.
116CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
On présente dans ce qui suit les modèles ARCH et GARCH ainsi que leurs principales
propriétés.
Une façon simple d’illustrer les dangers liés à l’utilisation des rendements au carré pour
mesurer la volatilité consiste à simuler un AR(1) et à représenter le carré du processus,
en le confrontant à une série financière. L’illusion est parfaite : on aurait presque l’im-
pression que les deux séries sont issues du même processus. C’est ce qu’on observe sur
la figure 5.23.
Enfin, on est également à la recherche de modèles statistiques pour les rendements qui
permettent de générer de la leptokurticité, i.e. des distributions à queues épaisses. On
remarque en général que les séries financières présentent une kurtosis supérieure à 3,
signifiant simplement que les queues de distribution de ces séries sont en général plus
épaisses que celles de la loi normale.
4 e−04
carré
0 e+00
0 50 100 150
Temps
ARMA au carré
12
8
carré
4
0
0 50 100 150
Temps
épaisse que la queue haute). Les densités empiriques de ces trois lois sont présentées
en figure 5.24. La VaR calculée pour différents seuils permet de se faire une idée plus
précise des risques liés à calculer une VaR sur une hypothèse de gaussiannité des ren-
dements, quand ceux-ci admettent des queues épaisses et/ou asymétriques.
La table 5.2 fournit une autre illustration de ces propriétés des séries financières : on
présente espérance, écart type (volatilité), skewness et kurtosis des rendements de l’in-
dice DAX ainsi que de divers call européens de strike différents, ayant tous le DAX
pour sous-jacent. On remarque tout à fait l’existence de queues épaisses ainsi que d’une
asymétrie.
Notons pour terminer cette section fourre-tout que ce qui vient d’être dit n’a pas que des
implications pour le risk-management, mais aussi pour l’asset pricing. A peu de choses
pret, comme le souligne très justement Cochrane (2002), le prix d’un actif est toujours
et partout une espérance sous loi risque neutre actualisée. Dans le cas où les rendements
sont gaussiens, la formule de Black Scholes semble tenir et est largement utilisée dans
de nombreux domaines, dont les options sur action et les produits de taux. Cependant,
l’existence des faits styliés qui viennent d’être mis en avant peut conduire à d’impor-
tantes erreurs de pricing : le modèle BS avec rendements GARCH introduit par Duan
et par Heston et Nandi semble conduire à des prix d’option plus proche de la réalité que
les prix BS standards. Ceci a d’importantes implications en terme d’asset management.
Densités
0.5
0.4
0.3
Density
0.2
0.1
0.0
−15 −10 −5 0 5 10
Fig. 5.24 – Densité de la loi de Laplace (rouge), normale (bleu) et Student (vert)
Il s’agit donc encore de processus applicables à des séries préalablement centrées, comme
dans le cas des ARMA. Notons que les séries des rendements sont théoriquement na-
turellement centrées : il s’agit simplement d’une conséquence de la martingalité des prix.
La variance conditionnelle n’a plus rien à voir avec celle des ARMA :
p
V[xt |ht ] = V[ ht t |ht ] (5.193)
= ht V[t |ht ] (5.194)
= ht E[2t |ht ] (5.195)
| {z }
=1 par hypothèse
= ht (5.196)
Ceci tient simplement au fait que xt soit naturellement un processus centré. Si ces
modèles semblent d’un abord pratiques, il n’en reste pas moins qu’il produisent natu-
rellement des erreurs de mesure sur la volatilité.
Pour ce qui de la variance, il est possible de procéder par récurrence. Utilisons tout
d’abord la loi de la décomposition de la variance :
Ultime propriété d’un processus ARCH(1), il est possible de montrer que le carré du
processus admet une représentation AR(1). On suit ici ce qui en est dit dans Poon
(2005)[Chapitre 4]. Notons νt la différence entre x2t et ht . On a alors :
ht = ω0 + ω1 x2t−1 (5.218)
⇔x2t − νt = ω0 + ω1 x2t−1 (5.219)
⇔x2t = ω0 + ω1 x2t−1 + νt (5.220)
On retrouve ainsi un processus AR(1) sur les carrés des résidus. Ceci a plusieurs impli-
cations pratiques :
5.3. LES MODÈLES ARCH-GARCH 121
– D’une part, un processus ARCH(1) ne semble pas saisir de façon adéquate les proces-
sus de volatilité financière : on a vu lors des applications des processus ARMA que
la volatilité de certains actifs semble présenter une structure plus proche des ARMA
(retour à la moyenne en cas de choc importants) que des AR. Il sera donc nécessaire
de complexifier légèrement la chose, afin d’accomoder cette caractéristique empirique.
E[xt ] = 0 (5.225)
Pour ce qui est la variance, il est ici nécessaire de déterminer, comme précédement,
une formule de récurrence pour parvenir finalement à exprimer la variance en passant
àP
la limite. On ne refait pas ici les calculs : on se contente de fournir le résultat. Si
| pi=1 ω1,i | < 1, alors la variance du processus existe et est de la forme :
ω
V[xt ] = Pp0 (5.226)
1− i=1 ω1,i
séries.
Il suffit donc de prouver que la kurtosis d’un processus ARCH est supérieure à 3. La
preuve est simplissime et s’appuie sur le lemme de Jensen.
La preuve d’établit comme suit, en se souvenant que f (x) = x2 est bien convexe : soit
xt un processus ARCH(p) donné. Sa kurtosis est alors :
E[x4t ]
Ku (x) = (5.228)
E[x2t ]2
Cette dernière égalité tient au fait que t suit une N (0, 1) donc la kurtosis est égale à 3.
Dans la mesure où sa variance est égale à 1, constatons simplement que le dénominateur
de la kurtosis est toujours égal à 1 pour une N (0, 1). D’où le résultat proposé. En
utilisant Jensen, on a alors :
E[x4t ] E [ht ]2
Ku (x) = > 3 =3 (5.239)
E[x2t ]2 E [ht ]2
CQFD.
5.3. LES MODÈLES ARCH-GARCH 123
Ceci est bien évidement vrai dans le cas où les séries sont centrées, comme c’est le
cas pour les processus ARCH. Pour montrer que les processus GARCH ne sont pas
asymétriques, il suffit de montrer que le numérateur est nul dans le cas d’un ARCH(p) :
Non, les processus ARCH ne permettent pas de prendre en compte tous les faits stylisés
de la finance : seule la leptokurticité est prise en compte.
où t ∼ N (0, 1). Comme précédement, on donne les moments conditionnels, à partir des-
quels on fournira finalement les moments non conditionnels par itération. L’espérance
conditionnelle du processus est la suivante :
où Ft est la filtration engendrée par les valeurs passées de xt , de x2t et de ht . La variance
conditionnelle est alors :
p
V[xt |Ft ] = V[ ht t |Ft ] (5.249)
p
= ht V[t |Ft ] (5.250)
= ht (5.251)
124CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Pour ce qui est de la variance non conditionnelle, on procède comme précédement, par
itération. On sait que :
Or, on sait que E[x2t−1 ] = E[ht−1 ]. On peut donc réécrire la précédente égalité de la
façon suivante :
Comme précédement, à la condition que |ω1 + ω2 | < 1, la série engendrée admet une
limite (la variance existe) conduit à :
ω0
E[x2t ] = (5.266)
1 − (ω1 + ω2 )
On retrouve les mêmes propriétés de leptokurticité que pour les processus ARCH. On
le montre de la même façon que précédement, en utilisant le lemme de Jensen. En
5.3. LES MODÈLES ARCH-GARCH 125
Si les GARCH sont à même de générer de la leptokurticité, ils sont en revanche incapable
de générer de l’asymétrie. On le montre de la même façon que précédement, en utilisant
la loi des espérances itérées :
p
E[x3t ] = E[( ht t )3 ] (5.273)
p 3
= E[ ht E[3t |ht ]] (5.274)
p 3
= E[ ht × 0] (5.275)
=0 (5.276)
avec t ∼ N (0, 1). Cette généralisation des GARCH(1,1) ne sera pas utilisée par la
suite : le lecteur soucieux d’aller jusqu’au bout des calculs pourra appliquer ce qui a
été fait précédement pour les ARCH(p).
Sur le plan technique, l’estimation des ARCH n’est guère plus complexe que les ARMA :
là encore, le processus xt n’est pas i.i.d. et il est nécessaire de travailler conditionnel-
lement au passé de xt et de ht , afin d’obtenir un processus i.i.d.. On connait les mo-
ments conditionnels d’un processus ARCH : celui-ci est conditionnellement gaussien,
d’espérance nulle et de variance égale à ht . Le problème ici, comme dans le cas des
MA, est que la variance est ”inobservable” : on a besoin des paramètres du proces-
sus ARCH pour obtenir la variance conditionnelle du processus. Comme précédement,
la seule façon d’estimer ce type de processus est d’écrire la log-vraisemblance, puis
de déterminer l’ensemble des dérivées dont on a besoin pour l’implémentation des
méthodes numériques décrites au chapitre 5.
x2
En remarquant dans le précédent réarrangement que htt = 2t , on trouve alors qu’à
l’optimum, le score est nul. On rappelle que le score est l’espérance de la dérivée de la
log-vraisemblance. En effet, on a, d’après ce qui vient d’être dit :
2
x
E t = E[2t ] = 1 (5.284)
ht
Afin d’obtenir les dérivées de la log-vraisemblance, il ne reste qu’à obtenir les dérivées
de la variance conditionnelle par rapport à chacun des paramètres. Elles se calculent
aisément :
∂ht
=1 (5.285)
∂ω0
∂ht
= x2t−1 (5.286)
∂ω1
5.3. LES MODÈLES ARCH-GARCH 127
On est alors en mesure de mettre en oeuvre ce qui a été présenté dans le cadre du cha-
pitre 5 : l’estimation se fait généralement par Newton-Raphson, en utilisant la matrice
BHHH comme approximation de la variance de l’estimateur des paramètres. Il est ce-
pendant nécessaire d’en dire un peu plus : on a fait l’hypothèse que l’erreur du modèle
(t )suivait une loi normale. Qu’en est il si tel n’est pas le cas ? Greene (2002)[Cha-
pitre 11] revient sur les implications de cette hypothèse : comme on l’a déjà précisé
précédement, l’estimation par maximum de vraisemblance conditionel avec bruit gaus-
sien conduit à l’obtention d’estimateurs consistent (i.e. qui converge presque surement
vers les vraies valeurs des paramètres). Dans ce cas, on parle de Pseudo Maximum de
Vraisemblance. Gourieroux et al. (1984) remarquent cependant que l’estimation de la
matrice de variance/covariance des estimateurs par l’inverse de la matrice d’information
de Fisher n’est pas bonne. Il est nécessaire d’uiliser l’approximation suivante :
V[θ] = H −1 F H −1 (5.289)
avec :
∂ 2 lnL
H = −E (5.290)
∂θ∂θT
∂lnL ∂lnL
F =E (5.291)
∂θ ∂θT
On remarque que la matrice 2 × 2 F est l’estimateur BHHH présenté au chapitre 5.
La matrice F est en général très simple à calculer, puisque prise en espérance. Il suffit
pour cela de déterminer les dérivées secondes de la log-vraisemblance. Dans le cas d’un
ARCH(1), les calculs conduisent généralement aux résultats suivants :
n
∂ 2 lnL
2
1X 1 xt
=− 2 −1 (5.292)
∂ω02 2
i=1 t
h2 ht
n
∂ 2 lnL 1 X x2t−1 x2t
=− 2 −1 (5.293)
∂ω0 ∂ω1 2
i=1
h2t ht
n
∂ 2 lnL 1 X x4t−1 x2t
=− 2 −1 (5.294)
∂∂ω12 2
i=1
h2t ht
On ne calcule que l’une des dérivées croisées, dans la mesure où le lemme de Monge
s’applique sans problème. En prenant l’espérance de ces dernières dérivées, et en re-
x2
marquant encore une fois que E[ htt ] = 1, on obtient la matrice F :
x2
− 21 ni=1 h12 − 21 ni=1 ht2
P P
F = 1 Pn x2t t
x4 (5.295)
− 2 i=1 ht2 − 21 ni=1 ht2
P
t t
128CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
On sait de plus calculer la matrice BHHH à partir des dérivées premières : on est alors
en mesure d’estimer un modèle ARCH sans trop de difficulté. Dans la réalité (et dans
R), l’estimation ne passe pas par Newton Raphson, mais par une méthode plus complexe
garantissant de meilleurs résultats : la méthode BFGS. Celle-ci n’est pas développée ici.
Il s’agit donc du même point de départ qu’on modèle ARCH(1). L’estimation d’un
modèle GARCH(1,1) nécessite cependant le calcul d’une dérivée première supplémentaire :
∂ht
= ht−1 (5.297)
∂ω2
On obtient alors les dérivées de la vraisemblance suivantes :
n
x2t
∂lnL 1X 1
=− 1− (5.298)
∂ω0 2 ht ht
i=1
n
1 X x2t−1 x2
∂lnL
=− 1− t (5.299)
∂ω1 2 ht ht
i=1
n
x2
∂lnL 1 X ht−1
=− 1− t (5.300)
∂ω2 2 ht ht
i=1
On laisse à titre d’exercice le calcul des dérivées secondes. On renvoie à VonSachs and
VanBellegem (2002) pour plus de précisions au sujet de l’estimation des GARCH. On
propose le code R suivant, stricte application de ce qui vient d’être dit : il s’agit d’un
code permettant d’estimer un GARCH(1,1) par Newton Raphson, avec matrice BHHH.
[Link]<-function(theta,x){
G=matrix(1,3,1)
n=nrow(x)
check=theta
j=1;
while(sum(G^2)>0.0001){
# Computation of sigma2
sigma2=matrix(theta[1,1],n,1);
for (i in 2:n){sigma2[i,1]=theta[1,1]+theta[2,1]*sigma2[(i-1),1]+theta[3,1]*x[(i-1),1]^2}
comp=[Link]((x^2-sigma2)/sigma2^2);
BHHH=cbind(comp[2:n,1],comp[2:n,1]*sigma2[1:(n-1),1],comp[2:n,1]*x[1:(n-1),1]^2);
H=(t(BHHH)%*%BHHH);
G[1,1]=sum(BHHH[,1]);
G[2,1]=sum(BHHH[,2]);
5.3. LES MODÈLES ARCH-GARCH 129
G[3,1]=sum(BHHH[,3]);
cat(j,"\n");
check=cbind(check,theta+solve(H)%*%G);
theta=theta+solve(H)%*%G;
j=j+1
}
var=sqrt(diag(solve(H)));
test=theta/var
return(list(theta=theta,check=check, test=test))
}
On étudie naturellement les ACF et PACF des rendements élevés au carré sur la fi-
gure 5.26. On remarque ce qui a été dit plus haut : l’autocorrélation dans les rende-
ments au carré décroit lentement et la PACF admet quelques valeurs significativement
différentes de 0. Ceci suggère naturellement l’estimation d’un modèle GARCH avec un
ordre faible. On choisit ici délibérément un GARCH(1,1), en utilisant la fonction R
fournie précédement.
La table [Link] fournit les estimations d’un modèle GARCH(1,1) sur nos données. On
remarque l’ensemble des paramètres est significatif à 95%. La figure 5.27 présente les
résidus du modèle GARCH. On rappelle que les résidus d’un GARCH sont données
par √xht et non par l’écart entre le modèle et les résidus. On constate que ces résidus
t
ne présentent presque plus de phénomènes de cluster de volatilité. En revanche l’étude
de la PACF et de l’ACF du carré des résidus (figure 5.28) conclue naturellement à
l’insuffisance de l’ordre du GARCH choisit : il subsiste encore de la corrélation.
0.05
−0.05
Time
0.4
0.0
Lag
0.05
−0.05
Lag
0.4
0.0
0 5 10 15 20 25 30
Lag
0.05
−0.05
0 5 10 15 20 25 30
Lag
rend
● ● ●
0.05
● ●●
● ● ●
● ● ●
●● ● ●●
●●● ●● ●
●●●● ●● ● ● ● ● ● ● ●
●●
● ● ● ● ●● ●● ●●●●● ● ●●●●● ●
●● ●● ● ●● ● ●
● ●●●● ● ● ●●● ● ● ●●
●● ●● ●● ● ● ● ●● ● ●●●● ● ● ● ● ● ●● ● ●● ●● ●● ● ● ●●● ●
●
●●● ● ●● ●
Series
● ●● ● ● ●●●● ● ● ●
●● ● ● ● ●
●● ●
●
● ● ●●●
● ●
●
●
●
Index
Residuals
● ● ● ●
● ● ● ● ●● ● ● ●
● ●
● ● ● ● ● ●
● ● ● ● ●● ●● ● ●●● ● ●
●● ●● ● ● ●●● ●●
2
● ●● ● ●●● ● ● ● ● ●●● ● ● ● ● ●● ●
●●● ● ●● ● ●● ●● ● ●● ● ● ●● ●●● ●●
●●●●● ●●
● ●● ● ● ● ●
●●● ●● ●● ●● ● ●
● ●● ●●● ● ●
● ● ●●● ● ●●
●●●
●
●●● ●
●● ● ●
●●● ● ● ● ●● ● ●● ● ● ● ●● ●●● ●●●
●● ●●●● ●● ●
●●●●●●● ●● ● ●● ●● ●●●●● ● ●● ● ●●●●●●● ● ●●
●●●●●
● ●● ● ●●● ●● ●
●●
●●●●
●
●● ●●
● ●●●
●●● ●
● ● ●●
●● ●
● ● ●● ● ●●
● ● ●
● ●●
●●
●● ●● ●●
●● ● ●●● ● ● ● ●●●●
● ● ● ●● ● ●
●
● ●●
●●●
● ●●●● ●
● ●●
●●●
●
●●
●● ●●
●●
●●
●
●
●
● ●● ●●
● ● ●● ● ●
● ● ●
● ●●●●
● ●● ●●● ● ●●
●●●●● ●●●
●
● ● ●● ● ●● ●●
● ●●
●
●
● ●
● ● ● ● ●● ●
● ●●● ●
●●●
● ●●● ●
●●● ●● ●● ● ●●● ●
● ●●
●●●●● ●
●●● ●●● ●●
●● ●
●
●● ●●●● ●
●● ● ●● ● ●●●
● ● ●● ●●
●●● ● ●●
● ● ●●● ● ●●● ●●
●● ●●●●● ●
●● ●●●● ● ●●
●●●●●
● ● ● ● ●
●● ●
● ● ●
● ● ● ●●●
●●●
●● ● ● ● ●● ● ● ● ● ● ● ●●● ● ●● ● ●● ● ●●
● ● ●● ●● ●●● ● ●●●
Series
●
●●
●●
● ● ●●●
●
● ●
●●
●●
●
●●●● ●
●●●●
● ●
● ● ●●
● ● ● ●
●●●
●●●
●● ●
● ●●●●● ●●● ● ●● ●
● ●●
● ●●● ● ● ●●
●●● ● ● ●●
●●●
●●
●●● ●
●●
● ●●●●●●
● ● ●●●●●
● ●
●●
● ●●
●●●●●
●●●
● ●●●● ●
●●● ●●● ●● ●
●
● ● ●●
●●●●
●●
●●●● ●
●●●●●
●
●●
● ●
●●
●●
●●
●●●
●●
● ●
●● ●●
● ● ● ●●●● ● ● ●● ● ● ● ● ● ●● ● ●● ● ●
● ● ●●● ● ● ● ● ● ● ●●● ● ●
● ●
● ● ● ●● ● ●●● ●
●●
● ● ● ● ●●●
● ●
●
● ●
● ●●● ●
●● ●●● ●
●● ●●●
0
●●●● ● ●● ●●●● ●
●●● ●●● ● ●
●● ●●● ●
●● ● ●● ● ●
● ●● ●● ●
● ●●●●● ● ● ● ●●● ●● ●●●● ● ●● ● ● ●
● ●●●● ● ● ●●●●● ●●●●● ●● ●
●●●●● ●● ● ● ●
● ●
●●● ● ● ● ●●● ●
●
● ●●●●●●
●●●●
●
●●●
●●
●●●●
● ●
●●●●●●●
●●●
● ●●●● ● ●● ●●
●●●● ●●●●●
●●
●●●● ●●●● ●●●● ● ●● ●
●●●●●●
●● ●●●●●●●●●●●●
● ●●
●●●
●
●●●●●
● ●● ●
●
● ●● ●● ●●
● ●● ●● ●
●●
●
● ●
●●●●
●●●●●●● ● ●●● ●●
●
●●●●
●
●●●●
● ●● ●●●● ●
●●●
●●●●●●
●●●●●● ● ●●
●
●●
●●● ●
●
●●●●
●●●
● ● ●●●●●
●●●● ●
●●
● ●●●●●●●●●
●●●●● ●
●●●● ●●
●●●
●
● ●
●●●●
●● ●●● ● ●● ●●●
●●
●●●
●●●● ● ● ●●
●● ●●● ●● ● ● ●●●●●●
● ● ● ●●●● ●●●●●● ●●●●●● ●
● ●● ●●
● ● ● ●●●●
●●●●●●●●
● ●●●●●
●
● ●● ●● ● ●● ● ●● ●● ●●●●
●●●●
● ●● ●● ●●
●●●● ●●●●● ●●● ●● ●● ●●●●● ●●
●● ●●●● ● ● ●
●● ● ●●●●●
● ●● ●
●● ●● ●●● ●●● ●●
●●
● ●●
●
● ●● ●● ● ● ●● ●● ● ●●● ● ●●●● ● ● ●●
● ● ●● ●
● ●● ●●● ● ● ● ●● ● ●● ● ●●
●● ● ●● ●● ● ● ●●●
●● ●● ●
● ●● ● ● ● ●● ● ●
●
●
●
●
● ● ●● ●● ●● ●●● ●
● ●● ●
●●
●●●
● ●● ● ● ●●● ●●● ● ● ● ● ●●
● ●●●●●
● ● ● ● ● ● ●● ● ●
● ● ● ●
● ●● ● ●● ●● ●● ●●
● ● ● ● ● ● ● ●●● ● ● ● ●
● ●●● ● ● ●● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●
● ● ●
●
−4
● ●
Index
0.4
0.0
0 5 10 15 20 25 30
Lag
0.4
0.0
0 5 10 15 20 25 30
Lag
Fig. 5.28 – ACF et PACF des rendements du DAX et des résidus GARCH
132CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
0.020
0.015
0.010
0.005
Time
Le tableau suivant fourni les erreurs des versions standards et GARCH ad-hoc générées
par la formule de Black-Scholes.
Une fois un modèle GARCH calibré, il est possible de procéder à une analyse plus fine
de la volatilité d’un sous-jacent. La volatilité est principalement utilisée dans le cadre de
la Value at Risk et dans le cadre des modèles d’option. L’avantage d’un modèle GARCH
est qu’il permet de construire une prévision naı̈ve de la volatilité future du marché, sur
une période courte. On présente ici quelques questions relatives à la prévision de la
volatilité du rendement d’un actif.
avec t ∼ N (0, 1), il est aisé d’obtenir une prévision en date t de la variance condition-
nelle en date t + 1 :
70
60
50
40
Prix option
30
20
10
0
0 50 100 150
Index
Fig. 5.30 – Prix du call DAX (rouge), prix BS standard (vert) et prix GARCH BS
ad-hoc
0.009
0.007
0 50 100 150
Index
30
0 10
0 50 100 150
Index
Fig. 5.31 – Prix théorique contre volatilité conditionnelle pour une option hors de la
monnaie.
5.3. LES MODÈLES ARCH-GARCH 135
Comme on l’a montré lors du calcul da variance non conditionnelle, en procédant par
récursion, on trouve la formule générale suivante :
h−1
X
E[ht+h |ht ] = ω0 (ω1 + ω2 )i + (ω1 + ω2 )h ht ) (5.312)
i=0
(5.313)
Comme dans le cas des ARMA, la prévision s’écrase rapidement contre la variance non
conditionnelle. Ainsi l’horizon prédictif des modèles GARCH est limité : ces modèles
sont bien mieux adaptés pour fournir une mesure de la volatilité, qui, rappelons le, est
par essence inobservable.
[Link].1 La VaR
La VaR est le benchmark le plus utilisé quand il s’agit de juger de l’exposition maxi-
male au risque de marché produit par une position donnée sur le marché. Le risque de
marché est un terme générique permettant de décrire plusieurs situations possibles sur
le marché, selon le montant des engagements produits. Une position conduisant à in-
vestir dans un seul et unique actif n’induit pas la même exposition au risque de marché
qu’une position prise sur différents actifs simultanément. Dans le cas où l’on détient
un portefeuille d’actifs, et en dépis des mécanismes bien connus de la diversification, la
dépendance existant entre les différents actifs du portefeuille fait peser sur sa rentabilité
une menace particulière, souvent présenté sous le titre de risque de corrélation.
La VaR peut être simplement définie comme la perte maximale liée à une position de
marché, sur une fenetre de temps particulière et pour un niveau de probabilité donné :
en supposant que l’on connaisse la loi du rendement du portefeuille, on est en mesure
de donné le montant maximal de perte que l’on peut réaliser, dans le cadre d’une pro-
babilité égale à 95%. Cela signifie simplement que dans 95% des cas, la perte ne devrait
pas dépenser cette mesure de risque. Evidement, la perte ne se juge pas en terme de
rentabilité ! Il est donc nécessaire de passer de la perte en terme de rentabilité à celle
en terme de prix, en passant à l’exponentielle.
136CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
L’idée est donc : connaissant la loi de t rt+h , i.e. des rendements sur la période de temps
allant de t à t + h (la fenêtre de temps), on cherche la chutte maximale que peut
connaitre le rendement du portefeuille sur ce même intervalle, avec une probabilité
égale à 90% ou 95%. Plus formellement, ceci revient à chercher :
Une fois ce quantilé obtenu (sous R, la function quantile permet d’obtenir une esti-
mation non paramétrique de ce quantile), il est alors possible d’obtenir une estimation
de la VaR en terme de prix, i.e. en terme de perte nette :
On obtient ainsi la perte nette maximale obtenue pour une probabilité égale à 95%.
Ceci amène plusieurs commentaires :
– Cette probabilité est en réalité une probabilité conditionnelle : sachant αPt , le mon-
tant net de l’engagement à la date t dans l’actif étudié, quel serait la perte maximale
rencontrée dans 95% des cas pour un horizon h.
– Il s’agit d’une VaR dans le cas d’une position courte : on a ici acheté le titre, et
seule une baisse de son cours (un enchainement de rentabilité négatives) peut nous
conduire à perdre de l’argent. On opère le raisonnement inverse dans le cas où l’on est
placé dans le cadre d’une position longue (vendeur à découvert) dans le titre étudié.
La perte dans ce cas viendra naturellement d’une montée du cours de l’actif. Il est
alors nécessaire de calculer la hausse maximale que le titre peut enregistrer dans 95%
des cas.
5.3. LES MODÈLES ARCH-GARCH 137
On rappelle simplement que la fonction quantile pour une variable aléatoire x est la
fonction telle que :
Comme présenté dans Tsay (2002), le calcul de VaR est loin d’être aisé. Il nécessite de
nombreux inputs :
Nombre de ces items sont en général fixés par le régulateur (accords Bâle II par
exemple). Notons que le calcul de la VaR sert aux comptables pour déterminer le
montant de provisions à passer pour couvrir le risque de marché. Il s’agit d’un montant
placé dans les capitaux propres (Provisions pour risques et charges) afin de faire face
à un décrochage violent mais temporaire du marché. Inutile de rappeler qu’un retour-
nement puissant et durable du marché conduit à des effets systémiques d’une ampleur
telle qu’il est impossible d’y faire face avec ce simple montant en poche. Ainsi, d’un
point de vu pratique, on fixe en général de façon institutionnelle un montant de VaR
qui sera passé en provision pour l’année à venir. Il sert ensuite à piloter le montant et le
sens des engagements dans les actifs en portefeuille : on détermine la VaR de l’ensemble
des titres détenus, et on vérifie que celle-ci ne dépasse pas le montant conventionnel-
lement fixé. Dans le cas où la banque est engagée dans différents marchés, et possède
donc différents desks, il est alors nécessaire d’allouer une part de la VaR conventionnelle
à chacun de ces desks.
On reviendra plus loin sur la méthodologie proposées par RiskMetrics (méthode EWMA),
dans la mesure où elle corresond à un modèle GARCH ”raffiné” (modèle GARCH
intégré).
Avec CAC ∼ N (0, 1). Les estimations obtenues sont les suivantes :
A l’aide de ces résultats, on peut déterminer une VaR pour l’horizon souhaité. Pour
une la VaR à un jour, on cherche :
CAC
≤ V aR1CAC = 5%
P rt+1 (5.324)
q
CAC CAC CAC
⇔P ht+1 t+1 ≤ V aR1 = 5% (5.325)
V aR1CAC
⇔P CAC
t+1 ≤ q
= 5% (5.326)
CAC
ht+1
suivant une loi normale centrée réduite, on connait le quantile à 5% : il vaut -1,64.
On en déduit donc :
V aRCAC
⇔q 1 = −1, 64 (5.327)
hCAC
t+1
q
⇔V aR1CAC = −1, 64 hCAC
t+1 (5.328)
5.3. LES MODÈLES ARCH-GARCH 139
La fin des calculs s’appuie sur des données numériques : on a besoin du prix du CAC en
date t ainsi que la prévision de la volatilité à la date t + 1. On a procédé à ces quelques
calculs sous R. On présente le résultat des opérations : on a calculé pour toutes les dates
la VaR à un jour. Elle est présentée en figure 5.32.
1800
1600
1400
1200
Index
Dans le cas où l’on souhaite obtenir une VaR pour un horizon t + h quelconque, les
choses se complexifient légèrement. On cherche à définir une borne basse pour Pt+h .
On travaille en rendements, i.e. sur :
Pt+h
r̃t+h = ln (5.332)
Pt
Il est alors nécessaire de travailler sur la distribution conditionnelle prédictive, i.e. sur
140CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
Il suffit donc de calculer la somme des variances conditionnelles pour déterminer ensuite
une VaR forward pour un horizon h. On représente la VaR et les moments où celle-ci
5.3. LES MODÈLES ARCH-GARCH 141
2600
2400
2200
2000 VaR pour CAC
Prix
1800
1600
1400
1200
Index
Avant toute chose, il est important de remarqué que la VaR d’un portefeuille, i.e. d’une
somme d’actifs pondérés, n’est pas jamais la sommes pondérée des VaR des différents
acifs. Ce qui rend le calcul de VaR d’un portefeuille complexe est la dépendance exis-
tante entre les différents actifs composant le portefeuille. Dans le cadre d’un modèle
142CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
On se propose de déterminer la VaR de la somme de ces actifs, non par calcul direct,
mais par simulation. Connaissant ρ, on est en mesure de simuler des corrélés, puis
de simuler les processus r1 et r2 conditionnellement à l’information disponible en t.
Enfin, on est en mesure de retrouver la VaR en terme de perte nette, comme on l’a fait
précédement, en prenant l’exponentielle de chaque rendements.
Une fois ceci fait, on estime la corrélation entre les résidus, en supposant que celle-ci
n’est différente de 0 qu’instantanément (pas de corrélation entre les résidus pour un
retard quelconque). On commence alors les simulations :
On a procédé à une application sur le portefeuille CAC+DAX pour une VaR à 10 jours.
On présente les résultats obtenus en figure 5.34. Les résultats semblent intéressants. Le
principal problème reposant sur le fait que lorsque l’on observe une chutte brutale du
prix, on est conduit à surestimer la VaR 10 jours après. Ce résultat est naturelle dans la
mesure où la VaR à 10 jours est, à un coefficient multiplicatif pret, le prix d’aujourd’hui.
L’intérêt principal d’une démarche basée sur des simulations est qu’elle permet, en
travaillant sur des séries univariées, d’incorporer la dépendance existant entre actifs.
Mieux, comme on le détaillera dans la section sur les GARCH, on est implicitement
amené à faire varier la covariance conditionnelle entre les deux actifs, en dépis du fait
que l’on ai fixé la corrélation/covariance entre les résidus une fois pour toute. Le prin-
cipal inconvénient est que, comme à chaque fois que l’on recourt à des simulations, le
temps de calcul est plus long que dans le cas univarié simple, présenté plus haut.
Notons pour terminer qu’il est possible d’effectuer l’ensemble de ce qui vient d’être dit
sur la base de calculs explicites. Dans la mesure où les deux actifs ont des rendements
conditionnellement gaussiens, et en constatant que :
V(X + Y ) = V(X) + V(Y ) + 2Cov(X, Y ) (5.348)
on est en mesure de décrire la loi conditionnelle de l’ensemble des deux rendements. En
effet, on sait que :
E[r1,t + r2,t |ht ] = 0 (5.349)
p p
V[r1,t + r2,t |ht ] = h1,t + h2,t + 2 h1,t h2,t ρ (5.350)
(5.351)
144CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
3800
3600
Portefeuille
3400
3200
3000
Index
Il est alors possible d’écrire la VaR pour l’horizon souhaité, sans trop de problèmes
(théoriques). Il évident que d’un point de vue pratique, l’utilisation de la corrélation
comme mesure de dépendance est plus qu’un pari risqué ! On rappelle simplement que
la corrélation n’est qu’une mesure de l’existance d’un lien linéaire entre actifs. Ces
questions sont à mon sens - pour l’instant - bien trop avancées pour être intégrées dans
ce cours. Les problèmes de mesure de dépendance sont abordées dans Embrechts et al.
(1999) et Embrechts et al. (2001). L’une des réponses actuelles (quoique ceci commence
à dater) à ces question de corrélations passe par les copules, outil statistique permettant
de travailler de façon plus aisée avec des fonctions de répartition multidimensionelles.
Le lecteur intéressé par les applications en finance de ce type d’outils lira Cherubini
et al. (2005) avec intéret.
Il existe de nombreux processus GARCH, existant pour des raisons diverses : prise en
compte de la lente décroissance de l’autocorrélation du processus de volatilité, prise
en compte du lien risque/rentabilité, prise en compte de l’asymétrie observée dans les
rendements ou liens avec des processus continus. On présente ici trois extensions pos-
sibles des processus GARCH : les modèles GARCH intégrés, les modèles GARCH-M,
les modèles GARCH asymétriques ainsi que le modèle GARCH de Heston.
5.3. LES MODÈLES ARCH-GARCH 145
[Link] GARCH-M
Les modèles GARCH-M furent initialement introduits par Engle et al. (1987) : l’ambi-
tion de l’article est de présenter un modèle permettant de rétablir le lien entre rentabilité
et rendement, cher à la finance de marché. L’idée est donc de faire dépendre la renta-
bilité conditionnelle du risque conditionnel lui-même, de façon linéaire le plus souvent.
Le modèle GARCH-M(1,1) est le suivant :
p
xt = λht + ht t (5.352)
ht = ω0 + ω1 x2t−1 + ω2 ht−1 (5.353)
avec t ∼ N (0, 1). Seule la première équation est modifiée par rapport à un GARCH
classique : on ajoute un terme multiplicatif de la variance, λht , que l’on appelle prime
de risque. On donne les moments conditionnels :
Ajoutons que les GARCH-M sont naturellement leptokurtiques, tout comme les GARCH.
Qu’en est il cependant de l’asymétrie du processus ? Plutot que de se lancer dans des
calculs longs et complexes, on a simplement procéder à des simulations : on simule
des GARCH-M avec des paramètres égaux pour toutes les simulations et on calcule à
chaque fois la skewness de la simulation. La figure 5.35 présente la densité estimée par
noyau de l’estimateur de la skewness. Les simulations ont été effectuées à l’aide d’un λ
négatif : on remarque que la série a de bonnes chances d’être asymétrique à gauche.
L’intéret de ces modèles est de permettre une mesure à chaque date de l’espérance
et de la variance des rendements. Naturellement, ceci fait penser aux modèles de type
146CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
2.5
2.0
1.5
Density
1.0
0.5
0.0
−8 −6 −4 −2 0
cor(1 , 2 ) = ρ (5.360)
On connait déjà les deux premiers moments conditionnels univariés, pour chacun des
deux actifs. Il ne reste plus qu’à déterminer la covariance conditionnelle entre r1 et r2 .
Il suffit de la calculer directement :
p p
Cov(r1,t , r2,t |h1,t , h2,t ) = Cov(λ1 h1,t + h1,t 1,t , λ2 h2,t + h2,t 2,t |h1,t , h2,t ) (5.361)
p p
= Cov( h1,t 1,t , h2,t 2,t |h1,t , h2,t ) (5.362)
p p
= h1,t h2,t Cov(1,t , 2,t |h1,t , h2,t ) (5.363)
p p
= h1,t h2,t ρ (5.364)
variance et covariance conditionnelle pour chaque titre. La corrélation entre les titres,
conditionnelle ou pas, reste la même :
1.0
0.5
0.5
Rentabilité
Rentabilité
● ● ● ● ● ● ● ● ● ●
● ●●
● ●● ● ●●
● ●●
0.0
0.0
−1.0 −0.5
−1.0 −0.5
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Risque Risque
1.0
0.5
0.5
Rentabilité
Rentabilité
● ● ● ● ● ● ● ● ● ●
● ●●
● ●● ● ●●
● ●●
0.0
0.0
−1.0 −0.5
−1.0 −0.5
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Risque Risque
Il est à noter ici que l’on a pas procédé à l’estimation de ces modèles sur des séries réelles,
et ce pour plusieurs raisons. L’une des principales est que le lien rentabilité/risque est
loin d’être stable selon les actifs étudiés. Il serait faux de croire que le lien espérance
variance est toujours et partout validé. Il existe différents effets microstructurels qui
peuvent expliquer ces phénomènes, au nombre desquels l’impact de la liquidité sur la
formation du prix du risque. Ce type de considération dépasse cependant largement
l’objet de ce cours.
Notons simplement que sur les données étudiées dans le cadre de la VaR, il est intéressant
de remarquer que le lien risque/rentabilité existe et semble robuste. On calcule la
corrélation entre les rendements et la variance conditionnelle telle qu’on l’obtient en
estimant un modèle GARCH(1,1). En notant ρ cette corrélation, la statistique de test
148CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
permettant de tester que ρ est bel et bien significativement différent de 0 est la suivante :
ρ √
t̂ρ = p n ∼ Tn−1 (5.366)
1 − ρ2
où Tn−1 est une distribution de Student à n − 1 degrés de liberté. En calculant le coef-
ficient de corrélation ainsi que la statistique de test présentée, on obtient les résultats
suivants :
CAC DAX
ρ 0,07229244 0,04553088
t̂ 3,125987 1,965682
Sur ces séries, il semble exister une corrélation significative et positive entre rendement
et volatilité. En reprenant les estimations obtenues au cours de l’étude sur la VaR, on
propose une estimation naive du paramètre λ pour les rendements du CAC et du DAX :
on estime par MCO l’équation de la moyenne du GARCH-M, à équation de la variance
connue. La table suivant fournit les résultats de l’estimation :
rendement (pas très loin de la différence d’ordre 1) est lui stationnaire (bruit blanc en
général).
ce processus n’est stationnaire (sa variance n’explose pas) que si ω1 + ω2 < 1. Dans
de nombreuses estimations de GARCH, on observe en général que la somme de ω1 et
ω2 très proche de 1, sans toutefois lui être supérieur. La table 5.4 présente l’estimation
de modèles GARCH(1,1) sur les différents indices européens contenus dans la base de
donnée EuStockMarkets disponible sous R. On remarque que la somme des coefficients
est très souvent proche de 1 : on parle souvent de processus intégré pour la volatilité.
Ceci fait au passage remarqué que la volatilité d’une actif financier n’est ni constante,
ni une variable lisse : il s’agit d’un processus agité, sujet à des changements de régime
le plus souvent. Il existe des modèles à changement de régime en séries temporelles, qui
ne sont pas développés dans le cadre de ce maigre cours. Un lecteur intéressé trouvera
une breve introduction dans Wang (2003)[Chapitre 5].
Ce qui suit s’inspire largement de Poon (2005)[Chapitre 4]. Un processus GARCH pour
lequel on a ω1 + ω2 = 1 est un processus non stationnaire en variance, dans la mesure
où sa variance non conditionnelle tend vers +∞. Le processus rt (les rendements) reste
cependant stationnaire au sens stricte. On dit que la volatilité conditionnelle suit un
processus IGARCH(1,1), et ni le moment d’ordre 2 ni le moments d’ordre 4 n’existent :
ils divergent tous vers +∞.
L’un des modèles de volatilité les plus utilisés en Risk Management est le modèle
EWMA de RiskMetrics. EWMA signifie Exponentially Weighted Moving Average : il
s’agit simplement d’un modèle spécifiant la volatilité comme une moyenne pondérée des
volatilités passées, avec des pondérations décroissant exponentiellement dans le temps.
D’une façon générale, le modèle s’écrit :
Pτ
β iσ2
σt+1 = i=0
2 Pτ t−i−1 i
(5.369)
i=0 β
τ est l’horizon de mémoire du modèle. Ce modèle, outre le fait qu’il soit utilisé dans de
nombreuses applications de RiskMetrics, a la particularité d’être proche d’un modèle
GARCH intégré. Pour prouver cette propriété, il suffit de faire une petit coup de pont
d’Avignon : on commence par retravailler l’expression d’un GARCH(1,1), puis on re-
vient sur le modèle EWMA pour étudier la similarité.
Ceci est vrai, même dans le cas où ω1 + ω2 = 1. L’important est que ω2 < 1, ce qui
semble vraisemblablement être le cas sur de nombreuses séries financières. On obtient
alors une forme intéressante de volatilité. Comparons ceci au modèles EWMA :
Pτ
β iσ2
σt+1 = i=0
2 Pτ t−i−1
i
(5.374)
i=0 β
Ici aussi, pour peu que β < 1, il est possible de passer à la limite pour trouver une
forme analytique à la variance :
∞
X
2
σt+1 = (1 − β) β i σt−i−1
2
(5.375)
i=0
Pour terminer, revenons sur le calcul de la VaR dans la perspective d’une volatilité
intégrée. Tsay (2002)[chapitre 7] propose une ecriture alternative au modèle RiskMe-
trics. La présentation diffère légèrement de celle qui en est fournie dans Poon (2005),
même reste globalement la même. Le modèle EWMA repose sur l’hypothèse que la
distribution conditionelle des rendements soit normale, d’espérance nulle et de variance
σt . La dynamique de la variance est décrite de la façon suivante :
σt2 = βσt−1
2 2
+ (1 − β)rt−1 (5.376)
où β < 1. On retrouve donc bien un modèle intégré, avec une variance non conditionnelle
qui n’est pas définie, dans la mesure où la somme des deux paramètres de la dynamique
de la variance est égale à 1 par construction. Il s’agit par conséquent d’un modèle
IGARCH(1,1), sans drift. Il est possible de fournir une prévision de la variance, comme
on l’a fait dans le cas d’un processus GARCH classique. On construit ces prévisions de
façon récursive, en rappelant que :
rt = σt t (5.377)
On sait que :
Dans un modèle intégré, la meilleure prévision de la volatilité que l’on puisse formuler,
sachant que l’information dont on dispose en t est réduite à la seule la volatilité d’au-
jourd’hui est σt . Une autre propriété remarquable découlant de ce qui vient d’être dit
est la suivante :
2 2 2
E[σt+h |σt ] = E[σt+h−1 + (1 − β)σt+h−1 (2t+h−1 − 1)|σt ] (5.383)
2
= E[σt+h−1 |σt ] (5.384)
En itérant la dernière équation, on trouve donc que la meilleure prévision que l’on
puisse formuler de la volatilité future pour un horizon t + h est la prévision que l’on
pourrait faire à un jour.
Dans le cas où la filtration ne se réduit pas à σt , mais peut être étendue à {σt , rt }, alors
la prévision à un jour diffère de la volatilité d’aujourd’hui :
2
E[σt+1 |σt , rt ] = βσt2 + (1 − β)rt2 (5.385)
152CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
On constate alors qu’un modèle intégré ne conduit pas à un retour de la prévision vers
le niveau moyen de l’échantillon : la volatilité prévue à long terme n’est ni plus moins
que celle d’aujourd’hui. Ceci est globalement en désaccord avec ce qu’on observe sur le
marché. En général, une phase de forte volatilité (bull market) est suivie d’une retour
au calme (bear market). Néanmoins, comme la table précédente le montrait, la volati-
lité, sans être complètement intégrée, n’est pas loin de l’être.
2 |σ , r ] = σ̂ 2 2 2 2
Fort de ces deux éléments (E[σt+h t t t+1 et σ̂t+1 = βσt + (1 − β)rt ), on est alors
en mesure de déterminer une VaR pour un horizon h quelconque, en reprenant P ce qui a
été dit au cours des sections précédentes. On cherche la loi conditionnelle de hi=1 rt+i ,
comme précédement. Sait faire de calcul, on sait qu’elle est gaussienne, d’espérance
nulle et de variance égale à la somme des variances, du fait de l’hypothèse d’absence
de corrélation sérielle dans les rendements. On a donc :
" h # h
X X
V rt+i |σt , rt = V [rt+i |σt , rt ] (5.386)
i=1 i=1
Xh
2
= σ̂t+1 (5.387)
i=1
2
= hσ̂t+1 (5.388)
(5.389)
On trouve naturellement :
√
V aRh = 1, 64 × hσt+1 (5.393)
On comprend alors mieux l’intéret de travailler sur des modèles permettant de générer
de l’asymétrie. Ils permettent de rendre compte d’un fait stylisé important en finance,
i.e. l’effet levier.
Il existe plusieurs modèles permettant de générer cet effet asymétrique. Ils reposent en
général sur l’introduction d’une variable indicatrice (i.e. valant soit 0, soit 1) qui prend
la valeur 1 dans le cas où les rendements sont négatifs. L’un des modèles les plus connus
est le modèle GARCH de Glosten, Jagannathan et Runkle (Glosten et al. (1993)), dit
GJR-GARCH. Il s’écrit de la façon suivante :
p
rt = ht t (5.402)
2 2
ht = ω0 + ω1 rt−1 + ω2 ht−1 + α1rt−1 <0 rt−1 (5.403)
154CHAPITRE 5. INTRODUCTION AUX MODÈLES DE SÉRIES TEMPORELLES
t ∼ N (0, 1). 1rt−1 <0 est une variable qui prend la valeur 1 lorsque les rendements en
date t − 1 sont négatifs. Dans ce cas là, on observe une volatilité ht à la date suivante
qui est plus importante que dans le cas où les rendements sont positifs, pourvu que
α >. On a encore un processus pour les rendements qui est conditionnellement centré.
La variance conditionnelle quant à elle dépend des rendements à la date passé :
2
E[ht |rt−1 > 0] = ω0 + ω1 rt−1 + ω2 ht−1 (5.404)
2
E[ht |rt−1 < 0] = ω0 + (ω1 + α)rt−1 + ω2 ht−1 (5.405)
L’une des utilisations que l’on on aurait envie de faire de ces processus GARCH est bien
évidement de les utiliser afin de valoriser des actifs financiers. Plus particulièrement, on
sait que les options sont des actifs dont le prix est particulièrement sensible aux varia-
tions de la volatilité. Un certain nombres d’articles de recherche se sont ainsi tournés
vers des méthodes permettant de valoriser des options vanilles (typiquement un call) en
utilisant des dynamiques GARCH pour les rendements du sous-jacent. Dans le cadre
de ce type de démarche, on butte sur deux types de difficultés :
Commençons par nous intéressé à la la limite possible d’un processus GARCH. Bien
évidement, il s’agit de la limite au sens financier du terme, i.e. pour un pas de temps
infinitésimal. Il ne faut surtout pas confondre la limite en finance, avec une limite
vers l’infini, souvent utilisée en statistique (par exemple dans l’énoncé de la loi des
grands nombres). Engle et Ishida (2002) montrent comment trouver les diffusions cor-
respondants à certains processus GARCH. On propose ici de développer le cas d’un
GARCH(1,1).
5.3. LES MODÈLES ARCH-GARCH 155
avec les hypothèses habituelles. On travaille ici sur la dynamique de la variance condi-
tionnelle, ht . Il s’agit d’une série d’astuces permettant de trouver par passage à la
limite le processus effectivement suivi par la variance conditionnelle. On commence par
réécrire la variance comme suit :
Cette dernière expression correspond à un cas particulier d’un processus général développé
par Engle et Ishida, nommé CEV-GARCH, par analogie avec les modèles CEV (Constant
Elasticity Variance) développés par la finance en temps continu. La version discrète d’un
CEV-GARCH est :
où Wt est un mouvement brownien standard. Dans notre cas, on trouve naturellement
comme limite d’un GARCH(1,1) la diffusion suivante :
La question suivante, et qui n’est pas abordée ici, est de parvenir à déterminer la dyna-
mique risque neutre des rendements. Heston[1993] propose le principe de Neutralisation
Locale au Risque. Un lecteur soucieux d’en savoir plus lira avec intéret le dernier cha-
pitre de Aboura (2005), ainsi que le papier de Heston. Ce point n’est pas développé
dans la mesure où il s’appuie sur des modèles à volatilité stochastique qui n’ont pas
encore été abordés.
Ce qui change cette fois-ci, c’est la façon d’écrire l’équation de la variance. Celle-ci est
écrite en terme de logarithme :
p
ln(ht ) = ω0 + ω1 ln(ht−1 ) + ω2 (|t−1 | − 2/π) − αt−1 (5.421)
p
avec t ∼ N (0, 1). Dans cep cas E[|t |] = 2/π, d’où le fait que l’on retire cette quantité
dans le terme ω2 (|t−1 | − 2/π) : on travaille sur une variable centrée. α est supposé
être positif (mais précédé d’un signe négatif), afin de permettre les effets levier tels
qu’on les a évoqué plus haut. On a ainsi :
p
E[ln(ht )|ht−1 , t−1 > 0] = (ω0 − ω2 2/π) + ω1 ln(ht−1 ) + (ω2 − α)t−1 (5.422)
p
E[ln(ht )|ht−1 , t−1 < 0] = (ω0 − ω2 2/π) + ω1 ln(ht−1 ) − (α + ω2 )t−1 (5.423)
5.3. LES MODÈLES ARCH-GARCH 157
Ainsi, selon les valeurs des paramètres, la réponse de la volatilité à des chocs positifs ou
négatifs peut différer, autorisant la modèlisation des fameux effets leviers. On se forge
rapidement l’intuition que lorsque t−1 < 0, on a alors une volatilité conditionnelle
plus importante que dans le cas contraire. On n’en dira pas davantage au sujet de ces
modèles, ce qui reste étant trop avancé.
rt = σt t (5.424)
t ∼ N (0, σ ) (5.425)
ht = lnσt2 (5.426)
A la différence des EGARCH, on propose une dynamique propre à la variance, avec une
loi propre. Wang (2003) propose simplement de donner à lnσt2 un pattern ARMA(p,q).
On a alors dans le cas d’un ARMA(1,1), la volatilité suivante :
Au final, on a donc deux bruits différents : un premier bruit lié à l’équation des ren-
dements eux-même, et un second bruit venant de la volatilité elle-même. Dans ce cas
précis, la loi de la volatilité est une loi log-normale, à support positif, tout comme le
χ2 que l’on obtient dans le cas d’un ARCH(1). On retrouve un processus proche d’un
ARCH(1) dans le cas où θ = 0.
Il s’agit d’un chapitre relativement léger : on fournit à chaque fois un certain nombre de
références pour tout lecteur soucieux d’approfondir ses connaissances (comme Pépino
par exemple).
Yt = m(Xt ) + t (6.1)
où t est un bruit blanc, d’espérance nulle. m(.) est une fonction arbitraire mais lisse
de xt . Plaçons nous en une date t, et supposons alors que X = x. Supposons également
que pour cette date, l’on dispose de n observations pour Yt . On peut alors écrire :
n n
1X 1X
yi = m(x) + i (6.2)
n n
i=1 i=1
159
160 CHAPITRE 6. BOITE À OUTILS STATISTIQUES
où ωt (x) est un poids qui est plus important pour les yt pour lesquels on a des xt proche
de x. L’estimation est ainsi déterminée à la fois par la distance entre xt et x et à la fois
par le poids accordé à cette distance.
Tout ceci peut sembler plutot abscond, et est principalement utilisé lors de l’estimation
des densités par noyau.
Kh (x − xt )
ωt (x) = 1 Pn (6.7)
n t=1 Kh (x − xt )
Cette paramétrisation permet évidement d’obtenir des poids dont la somme soit égale
à n. On obtient alors l’estimateur à noyau de Nadaraya-Watson, qui prend la forme
suivante :
Pn
Kh (x − xt )yt
m̂(x) = Pt=1n (6.8)
t=1 Kh (x − xt )
6.1. MÉTHODES NON-PARAMÉTRIQUES ET APPLICATION 161
Dans la pratique, il est nécessaire de spécifier une forme pour le noyau. On utilise dans
le cas le plus courant un noyau gaussien qui est défini par :
x2
1
Kh (x) = √ exp − 2 (6.9)
h 2π 2h
Tout l’art de l’estimation par noyau consiste alors à déterminer un h optimal. Tsay
(2002) présente quelques intuitions dans le cas d’un noyau d’Epanechnikov. On donne
ici simplement la règle que l’on se forge rapidement en utilisant ce type d’estimation :
un h trop faible conduit à ne pas assez lisser l’histogramme. Au contraire, un h trop
important lisse plus que de raison l’histogramme empirique.
[Link](x = x)
0.4
0.3
Density
0.2
0.1
0.0
−4 −2 0 2 4
Voici un exemple de fonction R permettant d’appliquer la méthode des noyaux. Elle est
relativement simple et nécessite d’optimiser le h à taton.
kernel<-function(X,n,h,min,max){
k=length(X);
support<-seq(min,max,length=n);
[Link]<-numeric(n);
for (i in 1:n){[Link][i]=1/(k*h)*sum(1/sqrt(2*pi)*exp(-((X-support[i])/h)^2))};
return(list([Link]=[Link], t=support))
}
162 CHAPITRE 6. BOITE À OUTILS STATISTIQUES
La figure 6.2 présente l’estimation de la densité d’une loi normale centrée réduite dans
le cas où h est mal choisi. On a utilisé le code précédent pour construire ces graphiques.
h=5 h=0.05
0.35
0.075
0.30
0.25
0.070
Loi normale centrée réduite
0.20
0.065
0.15
0.10
0.060
0.05
0.055
0.00
−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3
Support Support
On travaille à présent sur une matrice X de taille M(n, p), contenant les n observa-
tions des p séries que l’on souhaite étudier. Une ACP propose de déterminer un nombre
réduit de facteurs qui sont des combinaisons linéaires des éléments de la matrice X, de
façon à expliquer les liens existants entre les différentes séries. Ce faisant, on parvient
en quelques sortes à expliquer la matrice de corrélations de ces p séries. L’important
dans cette méthode est de permettre la construction d’un nombre de facteurs qui soit
6.2. ANALYSE DES DONNÉES 163
Y = X − 1n g T (6.12)
Il est alors aisé d’obtenir la matrice de variance/covariance des séries. Il s’agit simple-
ment de :
1 T
V = X X − gg T (6.13)
n
1
= Y TY (6.14)
n
La matrice des observations centrées réduites s’écrit donc naturellement :
1
0 ... 0
s1
0 1
s1 ... 0
Z = Y D1/s , avec D1/s = . (6.15)
.. .. ..
.. . . .
1
0 ... 0 sp
où D1/s est la racine de l’inverse de la première diagonale de la matrice V, i.e. une ma-
trice composée de l’inverse des écart-types de chaque variable. La matrice de corrélation
des variables s’obtient alors de façon aisée :
1 T
R= Z Z (6.16)
n
164 CHAPITRE 6. BOITE À OUTILS STATISTIQUES
Il est possible de procéder à une ACP sur la matrice V ou sur la matrice R. Cepen-
dant, l’utilisation de V présente le désavantage de ne pas fournir des facteurs stables à
toute combinaison linéraire de chaque variable. En clair : si on modifie de façon linéaire
une variable particulière, l’ACP sur V ne fournira pas ex-post les mêmes facteurs. On
préférera alors utiliser R pour l’ACP : R est insensible à toute transformation linéaire
variable par variable. Ceci revient en réalité à réaliser une ACP sur V, en utilisant une
métrique particulière, i.e. une mesure de distance entre observations, i.e. la métrique
définie par D1/s . On en dira pas plus sur ce point plus théorique que pratique.
On appelle intertie totale de X la moyenne pondérée des carrés des distances des
observations au centre de gravité. On a donc :
n
1 X T
Ig = Xi,. − g T D1/s Xi,. − g T (6.17)
n
i=1
n
1X T
= Zi,. Zi,. (6.18)
n
i=1
L’inertie peut être vue comme une sorte de variance totale de X : il s’agit de l’écart entre
chaque observation et le centre de gravité. Cette notion présente un certain nombre de
propriétés intéressantes. On commence par réécrire l’inertie de la façon suivante :
n p
1 XX 2
Ig = Zi,j (6.19)
n
i=1 j=1
p n
X X 1 2
= Zi,j (6.20)
n
j=1 i=1
p n
X 1 X (Xi,j − gj )2
= (6.21)
n
j=1
s2j i=1
p n
X 1 1X
= (Xi,j − gj )2 (6.22)
j=1
s2j n i=1
| {z }
s2j
p
X s2j
= (6.23)
j=1
s2j
= T r(R) (6.24)
=p (6.25)
Avec la métrique utilisée, l’inertie totale est toujours la même : elle est égale au nombre
de variables présentes dans X. L’idée sous-jacente à l’ACP est de proposer une méthode
permettant de résumer un grand nombre de variables en un nombre réduit de facteurs,
ces facteurs ayant été composés suivant une règle basée sur l’inertie. Le critère pour
former les facteurs est simple : on cherche k facteurs orthogonaux, de façon à former
une nouvelle base orthogonalisée permettant de représenter nos p variables. On cherche
ces facteurs de façon à ce que l’inertie de X dans cette nouvelle base soit la plus im-
portante possible.
6.2. ANALYSE DES DONNÉES 165
P2 = P (6.26)
P T D1/s2 = D1/s2 P (6.27)
fi = Zui (6.34)
où λi est la valeur propre associée au vecteur propre ui . On obtient ces derniers résultats
car par définition, une valeur propre est telle que :
Rui = λi ui (6.41)
166 CHAPITRE 6. BOITE À OUTILS STATISTIQUES
uT
i ui = 1 (6.42)
Terminons par quelques remarques : on a vu tout d’abord que la variance d’un facteur
est égale à la valeur propre associée au vecteur propre permettant d’obtenir ce facteur.
Ces facteur étant orthogonaux, la variance de la somme des facteurs est la somme de
la variance des facteurs, autrement dit la somme des valeurs propres. On a donc :
k k
!
X X
V fi = V (fi ) (6.43)
i=1 i=1
Xk
= λi (6.44)
i=1
Un autre résultat classique de l’algèbre linéaire, la trace d’une matrice carré est égale
à la somme de ses valeurs propres. Dans notre cas, on a donc :
k
X
T r(R) = λi (6.45)
i=1
On trouve donc que l’inertie associée à X (du moins à sa version centrée réduite) est
exactement égale à l’inertie dans la nouvelle base. La différence tient au fait que dans
la nouvelle base, les facteurs sont ordonnés par part d’inertie décroissante : du fait
de l’alogorithme de maximisation permettant de trouver la solution à notre problème
(alogrithme qui n’a pas été présenté ici), les vecteurs propres que l’on détermine arrive
par valeurs propres décroissantes. Autrement dit, le facteur 1 est associé au vecteur
propre pour lequel on a la plus grande valeur propre.
Pour résumer, l’ACP revient à remplacer les variables composant X qui sont corrélées,
par de nouvelles variables, les composantes principales qui sont combinaisons linéaires
des variables composant X, non corrélées entre elles, de variance maximale et les plus
liées en un certain sens aux variables composant X. L’ACP est ce qu’on appelle une
méthode factorielle linéaire.
4 an 0,27 0,89 0,93 1,00 0,95 0,91 0,91 0,89 0,88 0,88 0,85 0,83 0,79
5 an 0,27 0,86 0,90 0,95 1,00 0,94 0,94 0,92 0,91 0,90 0,87 0,85 0,81
6 an 0,23 0,83 0,88 0,91 0,94 1,00 0,96 0,97 0,97 0,90 0,88 0,86 0,83
7 an 0,25 0,81 0,86 0,91 0,94 0,96 1,00 0,95 0,95 0,91 0,86 0,84 0,82
8 an 0,25 0,79 0,84 0,89 0,92 0,97 0,95 1,00 0,98 0,90 0,89 0,87 0,85
9 an 0,26 0,77 0,83 0,88 0,91 0,97 0,95 0,98 1,00 0,91 0,90 0,88 0,86
10 an 0,23 0,80 0,85 0,88 0,90 0,90 0,91 0,90 0,91 1,00 0,90 0,88 0,86
15 ans 0,24 0,77 0,81 0,85 0,87 0,88 0,86 0,89 0,90 0,90 1,00 0,98 0,94
20 ans 0,22 0,74 0,79 0,83 0,85 0,86 0,84 0,87 0,88 0,88 0,98 1,00 0,94
30 ans 0,24 0,70 0,75 0,79 0,81 0,83 0,82 0,85 0,86 0,86 0,94 0,94 1,00
167
168 CHAPITRE 6. BOITE À OUTILS STATISTIQUES
Les deux graphiques suivantes présentent les valeurs propres et les vecteurs propres
extraits de la base de données. Idéalement, on ne souhaite retenir que quelques uns de
ces facteurs pour décrire le comportement de la courbe des taux. La méthode usuelles
est la méthode du coude : on ne retient que les valeurs propres précédant la formation
d’un coude sur le graphique. Ici, on ne retiendrait que 3 ou 4 facteurs. On considère
en général que la courbe des taux est guidée par trois facteurs, dont le premier étant
la source principale de ses mouvements. Ce facteur est en général identifié comme la
politique monétaire.
Valeurs propres
●
10
8
6
valeur
4
2
●
●
●
● ● ● ● ● ● ● ● ●
0
0 5 10 15 20 25 30
maturities
Factor 1
Factor 2
●
Factor 3
0.5
●
Factor 4
●
Correlations
● ●
0.0
●
● ●
●
●
●
−0.5
●
−1.0
0 5 10 15 20 25 30
Maturities
Berndt, E. K., Hall, R. E., Hall, B., and Hausman, J. A. (1974). Estimation and infe-
rence in nonlinear structural models. Annals of Economic and Social Measurement,
3 :653–665.
Black, F. and Scholes, M. (1973). The Pricing of Options and Corporate Liabilities.
Journal of Political Economy, (81) :637–654.
Bollen, B. and Inder, B. (2002). Estimating Daily Volatility in Financial Market Utili-
zing Intra Day Data. Journal of Empirical Finance, 9 :551–562.
Cherubini, U., Luciano, E., and Vecchiato, W. (2005). Copula Methods in Finance.
Wiley.
Embrechts, P., Lindskog, F., and McNeil, A. (2001). Modelling dependance with copulas
and applications to risk management. Preprint ETZH.
Embrechts, P., McNeil, A., and Straumann, D. (1999). Correlation and dependency
in risk management : properties and pitfalls. Departement of Mathematik , ETHZ,
Zürich, Working Paper.
169
170 BIBLIOGRAPHIE
Engle, R., Lilien, D., and Robbins, R. (1987). Estimating Time Varying Risk Premia
in the Term Structure : the ARCH-M Model. Econometrica, 55 :391–407.
Glosten, L., Jagannathan, R., and Runkle, D. (1993). On the Relation between the
Expected Value and the Volatility of the Nominal Excess Return on Stocks. Journal
of Finance, 48 :1779–1801.
Gourieroux, C., Montfort, A., and Trognon, A. (1984). Pseudo Maximum Likelihood
Methods : Applications to Poisson Models. Econometrica, 52 :701–720.
Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, M., and Teller, E. (1953). Equa-
tions of state calculations by fast computing machines. Journal of Chemical Physics,
21 :1087–1092.
Münk, C. (2004). Asset Pricing Theory. Lectures Notes for PhD Students.
Parkinson, M. (1980). The Extreme Value Method for Estimating the Variance of the
Rate of Return. Journal of Business, 53 :61–65.
Poon, S.-H. (2005). A Practical Guide to Forecasting Financial Market Volatility. Wiley
Finance.