Econométrie R
Econométrie R
Econométrie R
net/publication/371510150
CITATIONS READS
2 3,063
1 author:
Sami Mestiri
Faculté des Sciences Économiques et de Gestion de Mahdia
66 PUBLICATIONS 167 CITATIONS
SEE PROFILE
All content following this page was uploaded by Sami Mestiri on 05 April 2024.
Sami Mestiri
Faculté des sciences économiques et de gestion Mahdia
mestirisami2007 ∂ gmail.com
Table des matières
1
3 Les modèles à décalages temporels 32
3.1 Les modèles linéaires autorégressifs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.1.1 Formulation générale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.1.2 Test d’autocorrélation et méthodes d’estimation . . . . . . . . . . . . . . . . . 33
3.2 Les modèles à retards échelonnés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2.1 Formulation générale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2.2 Retard moyen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.3 Détermination du nombre de retards . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Distribution infinie des retards . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3.1 Modèle de Koyck (progression géométrique) . . . . . . . . . . . . . . . . . . . 38
3.3.2 Modèle de Solow (distribution de Pascal) . . . . . . . . . . . . . . . . . . . . 39
3.4 Exemples de modèles dynamiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4.1 Le modèle d’ajustement partiel . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4.2 Le modèle d’anticipation adaptative . . . . . . . . . . . . . . . . . . . . . . . . 41
4 La Cointégration 43
4.1 Principe de la cointégration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Test de cointégration : Test de Engel et Granger . . . . . . . . . . . . . . . . . . . . . 44
4.3 Le modèle à correction d’erreur (ECM) . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2
6.3.1 Le modèle à effets fixes individuels . . . . . . . . . . . . . . . . . . . . . . . . 57
6.3.2 Le modèle à effets aléatoires individuels . . . . . . . . . . . . . . . . . . . . . . 62
3
10 Les modèles de séries temporelles 107
10.1 Définition d’une série temporelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
10.1.1 L’opérateur retard-différence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
10.1.2 Fonction d’auto corrélation simple . . . . . . . . . . . . . . . . . . . . . . . . . 108
10.1.3 Fonction d’auto corrélation partielle . . . . . . . . . . . . . . . . . . . . . . . 108
10.1.4 Définition d’une série stationnaire . . . . . . . . . . . . . . . . . . . . . . . . . 109
10.1.5 Propriété d’un processus stationnaire . . . . . . . . . . . . . . . . . . . . . . . 109
10.1.6 Le problème de la non stationnarité . . . . . . . . . . . . . . . . . . . . . . . 110
10.1.7 Nature de la non stationnarité . . . . . . . . . . . . . . . . . . . . . . . . . . 110
10.1.8 Tests de racine unitaire : les tests de Dickey-Fuller . . . . . . . . . . . . . . . . 111
10.2 Modélisation d’une série temporelle stationnaire . . . . . . . . . . . . . . . . . . . . . 113
10.2.1 Représentation de Wold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
10.2.2 Modèle autorégressif d’ordre p . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
10.2.3 Processus moyenne mobile d’ordre q . . . . . . . . . . . . . . . . . . . . . . . 118
10.2.4 Modèle auto régressif moyenne mobile ARMA (p,q) . . . . . . . . . . . . . . 120
10.2.5 Le processus autorégressif à moyenne mobile intégrée (ARIMA) (p,d,q) . . . . 121
10.3 La méthode de Box et Jenkins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
10.3.1 Transformation des données . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
10.3.2 Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
10.3.3 Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
10.3.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
4
12.4 Estimation du modèle à correction d’erreurs . . . . . . . . . . . . . . . . . . . . . . . 142
12.4.1 Exemple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
12.5 Généralisation à k variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
12.5.1 La cointégration entre K variables . . . . . . . . . . . . . . . . . . . . . . . . 144
12.5.2 Estimation du modèle à correction d’erreur . . . . . . . . . . . . . . . . . . . 145
12.6 Le modèle à correction d’erreur vectoriel . . . . . . . . . . . . . . . . . . . . . . . . . 146
12.6.1 Tests de relation de cointégration . . . . . . . . . . . . . . . . . . . . . . . . . 147
12.6.2 Test de la trace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
12.6.3 Test de la valeur propre maximale . . . . . . . . . . . . . . . . . . . . . . . . . 149
12.7 Synthèse de la procédure d’estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
12.7.1 Application avec R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
5
Avant propos
6
Chapitre 1
Dans les sciences sociales, et particulièrement en économie, les analystes étudient des phénomènes
afin de mieux comprendre la nature et le fonctionnement des systèmes économiques. L’objectif de
l’économétrie est de permettre aux agents économiques (ménages, entreprises, État) d’intervenir de
manière plus efficace.
On peut définir un modèle économétrique comme une présentation formalisée d’un phénomène
sous forme d’équations dont les variables sont des grandeurs économiques. L’objectif d’un modèle
est de représenter les traits les plus marquants d’une réalité qu’il cherche à styliser (interpréter en le
simplifiant).
L’économétrie est donc l’outil à la disposition de l’économiste qu’il lui permet de comprendre
et d’expliquer des phénomènes. Le modèle est l’outil que le modélisateur utilise lorsqu’il cherche à
infirmer ou de confirmer des théories qu’il construit. Pour ce faire, il émit des hypothèses et explicite
des relations.
7
Y(n,1) = X(n,k) β(k,1) + ε(n,1) , (1.2)
L’estimation de la spécification (1.1) nécessite préalablement la formulation d’un ensemble d’hy-
pothèses concernant le terme aléatoire (εi ) :
H 1 : E(εi ) = 0 les effets des facteurs autre que xi se compense.
2
H 2 : V (εi ) = σ Il s’agit d’hypothèse d’homoscèdasticité.
H 3 : Cov(εi , εj ) = 0 ∀i 6= j les erreurs sont indépendantes.
0
SCR = kY − Xβk2 = (Y − Xβ) (Y − Xβ). (1.3)
Le principe de la méthode des moindres carrés ordinaires consiste à minimiser la somme des
carrées résiduelles, par rapport aux paramètres inconnus du modèle. La résolution de ce programme
d’optimisation permet d’obtenir l’estimateur suivant :
0 0
β̂ = (X X)−1 X Y. (1.4)
D’après le théorème de Gauss-Markov, β̂ est le meilleur estimateur linéaire sans biais (BLUE, pour
Best Linear Unbiased Estimator). En effet, ce théorème prouve que β̂ possède la variance minimale
dans la classe des estimateurs linéaires centrés. La matrice de variance-covariance de β̂ a la forme
suivante :
0
V (β̂) = σ 2 (X X)−1 . (1.5)
Concernant l’estimation du paramètre de la partie aléatoire σ 2 , un estimateur centré de σ 2 est
calculé selon cette relation :
SCR
σ̂ 2 = , (1.6)
n−k
avec p est le nombre de paramètres à estimer.
8
ε ∼ N (0, σ 2 In ). (1.7)
Comme conséquence de la normalité des erreurs, nous obtenons :
SCR
∼ χ2 (n − k). (1.9)
σ2
Ceci permet d’effectuer les tests individuels à partir des statistiques distribuées selon la loi de Student
où la variance σ 2 est remplacée par son estimateur sans biais σ̂ 2 , il s’en suit que :
β̂ − βp
qp ∼ t(n − k). (1.10)
V (β̂k )
Dans ce qui suit, nous présentons deux types de tests individuels :
Le test d’a priori théorique comme par exemple β est égal à une constante. Les coefficients testés
prennent des valeurs réelles particulières selon l’appréhension
( théorique du modèle étudié. Dans ce
H0 : βk = β0
cas, nous testons les corps d’hypothèses suivants :
H1 : βk 6= β0
(n−p) β̂k −β0
On rejette H0 lorsque |tc | > t1− α avec tc = √ est l’estimation de la statistique du test à
2 V (β̂k )
(n−k)
partir de l’échantillon utilisé ett1− α
étant la valeur tablée.
2
Le test de Student des coefficients est généralement utilisé pour vérifier les hypothèses théoriques.
Mais, ce test n’est valable que lorsque l’hypothèse H2 est vérifiée c’est-à-dire le terme d’erreur est
vraiment un bruit blanc (ni autocorrelation ni hétéroscédasticité).
9
Nous introduisons ici les outils de modèle linéaire sous R, qui seront appliqués au cas de la régression
multiple. Les modèles linéaires sont des objets de classe lm que l’on crée par la fonction lm( ). Cette
fonction prend en argument une formule, puis éventuellement quelques options.
Le point important est la spécification du modèle par une formule. Pour le cas de la régression
linéaire simple, c’est assez simple : la régression de y par x s’écrit tout simplement :
lm(y ∼ x)
ou y et x sont des variables, vecteurs pour la régression, mais nous verrons plus tard que x peut
aussi être un facteur. Il est à noter que par défaut, les modèles incluent toujours l’intersection avec
l’origine.
Si on veut effectuer une régression de type y = a.x et non y = a.x + b, il faut explicitement enlever
l’intersection dans la formule :
lm(y ∼ x − 1)
Les variables à utiliser sont souvent dans un objet dataframe. La fonction lm peut prendre plusieurs
arguments, notamment :
-weights un vecteur de poids pour l’ajustement selon les moindres carrés,
- na.action une fonction pour traiter les données manquantes (NA = Not Available). Attention, le
fonctionnement est ici quelque peu différent de la fonction cor, puisqu’on passe une fonction et non
une chaîne de caractères. Plusieurs fonctions sont prédéfinies et doivent être utilisées :
- "na.fail" (défaut, équivalent à "all.obs") renvoie une erreur si une donnée manquante est rencon-
trée,
- "na.omit" enlève toutes les observations ayant au moins une donnée manquante (équivalent à "pair-
wise.complete.obs" dans le cas de la corrélation),
-"na.exclude" est similaire à la précédente, mais marque "NA" dans les sorties de fonctions.
- subset permet de travailler sur un sous-jeu de données : prend en argument un vecteur (ex : c(1 :10)
travaille seulement sur les 10 premières lignes, c(-2,-4) enlève les lignes 2 et 4, etc.)
y=c(10,8,7,12,11,13,15)
x1=c(5,3,2,6,7,8,10)
x2=c(7,13,16,10,18,22,30)
mod1=lm(y ∼x1+x2)
summary(mod1)
anova(mod1)
10
-summary imprime un résumé du modèle, avec coefficients et tests de significativité.
- print imprime une courte description du modèle.
- deviance calcule la somme des carrés résiduelle.
- predict (prend en argument un deuxième dataframe avec les mêmes variables explicatives que celui
utilisé pour ajuster le modèle) permet de prédire des valeurs à partir du modèle.
- AIC calcule le critère informatif d’Akaike(égal à −2.ln(L) + 2K, où L est la vraisemblance du
modèle et K le nombre de paramètres) . Peut prendre plusieurs modèles en argument.
- abline ajoute la droite de régression à un graphe.
i 1 2 3 4 5 6 7 8 9 10
cs 0.1 1.1 1.7 0.9 20.7 10.1 38.2 2.3 4.8 45.8
km 5 20 25 15 300 140 550 25 60 650
Commandes permettant de préparer les données :
cs=c(0.1, 1.1, 1.7, 0.9, 20.7, 10.1, 38.2, 2.3, 4.8, 45.8)
km=c(5, 20, 25, 15, 300, 140, 550, 25, 60, 650)
# Création d’un jeu de données R (data.frame) à partir d’un fichier
texte
consom = read.table("consom.txt")
# Nomme les variables V1 et V2 en km et cs
names(consom) =c("km","cs")
# Rend visible les vecteurs constituant consom
attach(consom)
csi = b + akmi + ei
11
-Le coefficient de qualité R2 (carré du coefficient de corrélation)
cov(x; y)2 var(ŷ)
R2 = r2 (x; y) = =
var(x)var(y) var(y)
Plus ce coefficient est proche de 1, meilleur est le modèle.
Le σ̂ 2 P 2
2 êi
σ̂ =
dll
avec ddl = degré de liberté = nombre d’observations - nombre de paramètres inconnus
Plus le σ̂ 2 , mesurant l’importance des résidus, est faible, meilleur est le modèle.
-Les erreurs standards
Elles doivent être faibles car elles mesurent l’incertitude des paramètres inconnus.
-Les tests realises sur les paramètres
Les tests de Student réalise sur chaque coefficient doivent être significatifs pour déduire que la variable
explicative a significativement un effet.
-Vérification de la qualité du modèle
Afin de vérifier si le modèle est bon, on doit visualiser :
-Les résidus par rapport aux valeurs prédites
Il ne doit pas exister de liaisons entre les points.
-Les valeurs prédites par rapport aux valeurs observées
Les points doivent être ajustes sur la première bissectrice.
12
Pour afficher, les éléments contenant les résultats d’une analyse, on utilise la fonction names().
names(regs)
"coefficients" "residuals" "effects" "rank"
"fitted.values" "assign" "qr" "df.residual"
"xlevels" "call" "terms" "model"
names(summary(regs))
"call" "terms" "residuals" "coefficients"
"aliased" "sigma" "df" "r.squared"
"adj.r.squared" "fstatistic" "cov.unscaled"
Il est aussi possible d’extraire les résultats fournis afin de pouvoir les réutiliser.
summary(regs)["sigma"]
$sigma
[1] 0.3839064
- La droite de régression obtenue pour ce modèle est : cs = 0, 0574 + 0, 0699 km. Une augmen-
tation de vitesse d’un km induit une augmentation de consommation de 0.067 litre.
- Ce modèle a un R2 très proche de 1 (0,9995) c’est à dire que la variation de consommation d’essence
(cs) est expliquée à 99.95% par la distance parcourue (km), c’est donc un très bon ajustement des
données.
- L’estimation de la variance (σ̂ 2 = 0.14738 ) obtenue est faible donc tous les points sont très proches
de la droite de régression.
- Les erreurs standards sont très faibles (proches de 0), montrant qu’il y a peu d’incertitude pour
l’estimation des paramètres.
- La P-Value du test sur le coefficient correspondant au km étant inférieure à 5%, il est possible de
dire que les km ont significativement un effet sur la consommation.
Il est donc possible de prédire la consommation de la voiture avec précision en fonction de la distance
à parcourir. Cela veut dire que la voiture consomme sept litres aux cent kilomètres.
Il est aussi intéressant de réaliser ses propres graphiques d’analyse. On peut représenter el nuage de
points des observations et la droite de régression en faisant :
plot(km,cs, main="Nuage de points des observations et droite de régression")
abline(regs)
Le nuage de points des valeurs prédites en fonction des valeurs observées
13
plot(predict(regs),cs, main="Nuage de points des valeurs prédites en fonction des valeurs observées")
Il est aussi possible d’extraire les résultats fournis afin de pouvoir les réutiliser.
Les fonctions residuals() et predict() permettent d’afficher les résidus calcules et les valeurs ajustées
de la régression.
round(residuals(regs),3) # arrondi à la troisième décimale
-0.307 -0.355 -0.105 -0.206 -0.328 0.256 -0.304 0.495 0.548 0.306
round(predict(regs),3) # arrondi à la troisième décimale
0.407 1.455 1.805 1.106 21.028 9.844 38.504 1.805 4.252 45.494
14
- De plus, le graphique représentant les consommations prédites en fonction des consommations ob-
servées montrent que les points sont alignés avec la première bissectrice ce qui permet de valider
ce modèle. En effet, les prévisions effectuées sont bonnes, elles doivent normalement être égales aux
valeurs observées c’est à dire que les points du graphique doivent être alignés sur la droite d’équation
"cs = predict( regs)". Dans ce cas, il semble que tous les points soient très proches de cette droite.
- Pour cet ensemble de raisons, il est possible d’affirmer que ce modèle est un très bon modèle.
15
H0F : Tous les coefficents valent zero
H0JB : Les residus suivent une loi normale
H0DW : Les residus ne sont pas auto-correles d’ordre 1
H0BG : Presence d auto correlation d’ordre superieur
H0BP : Heteroscedasticite des residus
H0GQ : Homoscedasticite des residus
H0R : Forme lineaire correcte
(19838.62 − 11010.44)/2
Fc = ) = 13.63
11010.44/34
16
1.3.2 Test de causalité de Granger
Afin de tester si X est une cause pour Y, il suffit de construire le modèle de la façon suivante :
k
X k
X
Yt = αj Yt−j + βj Xt−j + εt
j=1 j=1
(
H0 : β1 = β2 = ... = βk = 0
Puis effectuer le test d’hypothèse suivant :
H1 : β1 6= 0; β2 6= 0, ..., βk 6= 0
Pour faire ce test proprement, on fait appel à la statistique de Fisher habituelle sur un sous ensemble
de coefficients.
Enfin, il faut toujours tester la causalité dans l’autre sens, c’est-à-dire construire le modèle de la
façon suivante :Xt = kj=1 αj Xt−j + kj=1 βj Yt−j + εt
P P
Remarque : k qui correspond au nombre de périodes de retard est choisie par l’économètre en
fonction de la taille l’échantillon (en général 2, 4,8. . .).
Ou bien encore :
R2
c 1
F = 1−R2
n−k−1
Si F > Fαtabulé (1, n − k − 1) On rejette H0 c’est à dire le modèle est globalement significatif.
17
1.5 Prévision
Un modèle économétrique sert, en général, à prendre une décision dans la politique économique et à
établir des prévisions à moyen et long terme.
On distingue deux types de prévisions :
1. Prévision ponctuelle.
2. Prévision par intervalle
Soit le modèle vérifiant les hypothèses de GM : yi = β1 + β2 X2i + .... + βk Xki + ui i=1,. . ..,N
On veut prévoir yf /les valeurs prises par les exogènes dans le futur sont :(X2f , ..., Xkf ).
Exemple : Ci = a + bRi i=1980,. . .. . .,2006.
⇒b
c2007 = b
a + bR
c 2007 C’est la prévision ponctuelle de la consommation future.
q
α/2
Avec un risque α ∈ ]0, 1], on peut déduire un intervalle de confiance pour yf : yf = ybf ±tN −K V\
ar(ef )
avec ef est l’erreur de prévision :
b2 1 + Xf0 (X 0 X)−1 Xf
uf = yf − ybf Donc : V \
ar(uf ) = σ
18
Chapitre 2
Dans les sciences sociales, et particulièrement en économie, les analystes étudient des phénomènes
afin de mieux comprendre la nature et le fonctionnement des systèmes économiques. L’objectif de
l’économétrie est de permettre aux agents économiques (ménages, entreprises, État) d’intervenir de
manière plus efficace.
On peut définir un modèle économétrique comme une présentation formalisée d’un phénomène
sous forme d’équations dont les variables sont des grandeurs économiques. L’objectif d’un modèle
est de représenter les traits les plus marquants d’une réalité qu’il cherche à styliser (interpréter en le
simplifiant).
L’économétrie est donc l’outil à la disposition de l’économiste qu’il lui permet de comprendre
et d’ expliquer des phénomènes. Le modèle est l’outil que le modélisateur utilise lorsqu’il cherche à
infirmer ou de confirmer des théories qu’il construit. Pour ce faire, il émit des hypothèses et explicite
des relations.
L’objectif de ce chapitre est de présenter le modèle de la régression linéaire. Nous abordons tout
d’abord la notion du modèle statistique ainsi que nous envisageons une procédure d’estimation des
paramètres en étudiant les propriétés statistiques des estimateurs. La deuxième partie sera consacrée
aux problèmes particuliers liés au non respect des hypothèses.
Nous nous attachons particulièrement à deux formes classiques à savoir l’auto corrélation des
erreurs et l’hétéroscédasticité. L’étude de ces deux phénomènes nous permet de définir l’estimateur
des moindres carrés généralisés, utilisé lorsque la matrice de variances et covariances de l’erreur ne
réponds plus aux hypothèses classiques. Pour chacun de ces deux cas, nous présentons les méthodes
d’investigation : tests de détection, conséquences et procédures d’estimation.
19
Diagnostic : Une observation atypique se signale par un écart important sur le diagramme des
résidus. Lorsque des valeurs atypiques ont été trouvées, il faut commencer par vérifier les données
correspondantes, des erreurs ayant pu se glisser dans les chiffres. Si ce n’est pas le cas, il s’agit
en général d’observations particulières pour lesquelles on pouvait s’attendre à des singularités (par
exemple en octobre 1929, en mai 1968 ?.)
Solutions : On peut réestimer l’équation en supprimant la, ou les observations contestées. Et on
examine en particulier comment sont modifiées les différentes estimations. Cependant, une certaine
dispersion des résidus est la conséquence normale de la présence d’un aléa ; il serait absurde de vouloir
supprimer systématiquement toutes les observations un peu marginales !
2.2 Multicolinéarité
Il ne s’agit pas là d’un problème lié à l ?aléa, mais aux données disponibles pour l’estimation.
On suppose que sur les observations disponibles, X et Z sont exactement colinéaires (Z=3.X). La
0
matrice (X X) dans ce cas n’est plus inversible. Plus généralement, la multicolinéarité stricte entre
les variables explicatives rend la régression impossible (certains logiciels identifient et signalent cette
situation). Cette situation est néanmoins facilement évitable, car elle résulte le plus souvent d ?une
relation évidente entre explicatives envisagée, par exemple vouloir utiliser le budget civil, le budget
militaire et le budget total pour un ensemble de pays.
20
économétrique vraiment satisfaisante. Résultant d’un conflit entre la taille des données et celle du
modèle, il accepte deux types de solutions. On peut tenter de réduire le modèle, en éliminant une
variable mineure qu’on avait retenu dans le but de n’omettre aucun facteur explicatif, si du moins
elle peut être retirée sans altérer l’esprit du modèle. Plus de détail dans l’économétrie II.
21
a) Examen visuel des résidus
L’analyse graphique des résidus permet le plus souvent de détecter un processus de reproduction
des erreurs tel que :
-Les résidus sont pendant plusieurs périodes consécutives soit positifs, soit négatifs : auto corrélation
positif.
-Les résidus sont alternés : auto corrélation négatif.
Cependant le plus souvent, l’analyse graphique est délicate d’interprétation car le dessin de résidus
ne présente pas des caractéristiques toujours évidentes pour cela on doit effectuer un test statistique.
Le test de Durbin et Watson (DW) permet de détecter une autocorrélation des erreurs d ?ordre
1 selon la forme :
εi = ρεi−1 + ui
Le test d ?hypothèses est le suivant :H0 : ρ = 0 H1 6= 0. Le coefficient d’autocorrélation ρ est
significativement égal à 0.
La première étape consiste à estimer le modèle par la méthode MCO et à évaluer les résidus ei .
Ensuite à l’aide de ces résidus, on calcule la statistique de Durbin-Watson :
n
(ei − ei−1 )2
P
i=2
DW = n
e2i−1
P
i=1
22
Or −1 < ρ < 1 d’ou si ρ = −1 alors DW ' 4 et si ρ = 1 alors DW ' 0, enfin si ρ = 0 alors
DW ' 2.
Afin de tester l’hypothèse H0 , Durbin et Watson ont tabulé les valeurs critiques de DW au seuil
de 5% en fonction de la taille d’échantillon n et du nombre de variable explicative k. La lecture de
table permet de déterminer deux valeurs d1 et d2 compris entre 0 et 2 qui éliminent l’espace 0 et 4.
La règle de décision selon la valeur de la statistique (DW ) est représente dans le table (1)
Pour appliquer ce test, les conditions suivantes doivent être respectées : le modèle de régression
doit inclure une constante et les variables explicatives doivent être certaines. Le test de Durbin et
Watson est un test présomptif d ?indépendance des erreurs du fait qu ?il utilise les résidus ; de plus,
il ne teste qu ?une autocorrélation d ?ordre 1.
Exemple
Soit le modèle à trois variables explicatives :
Nous nous proposons de déceler une éventuelle autocorrélation d’ordre 1 des erreurs, pour cela
on demande :
1) d’estimer les coefficients du modèle ;
23
2) d’effectuer l’analyse graphique des résidus ;
3) de calculer la statistique de Durbin et Watson et d’effectuer le test ;
Solution : 1)Les résultats de l’estimation à partir de logiciel R sont les suivants :
(16)
Les trois coefficients des variables explicatives sont significativement différents de 0 (t0,05 = 2, 12),
(3,16)
le Fisher empirique est supérieur au Fisher lu (F0,05 = 3, 24), le modèle semble satisfaisant.
reg1=lm(Y∼X1+X2+X3)
summary(reg1)
2) L’analyse de ce graphique révèle des résidus qui semblent cycliques, ceci est symptomatique
d’une autocorrélation positive des résidus.
resd=reg1$residuals)
plot(resd, type = "l", main ="Résidus de Y") abline(h=0)
3) Les conditions d’utilisation du test de Durbin et Watson sont bien respectées : le modèle est
spécifié en série temporelle, le nombre d ?observations (n = 20) est supérieur à 15 et, enfin, le modèle
24
estimé comporte un terme constant. Le calcul de la statistique à partir des résidus est alors :
n
(ei − ei−1 )2
P
i=2 171.23
DW = n = = 1.053
162.49
e2i
P
i=1
Cette valeur est à comparer à celles lues dans la table de Durbin et Watson à n = 20 et k = 3 , soit
d1 = 1,00 et d2 = 1,68 . La valeur DW se situe dans la zone de doute (d1 < DW < d2 incertitude),
cependant à proximité immédiate de la zone de rejet de H0 , nous pouvons plutôt conclure à une
autocorrélation positive des résidus, donc à une présomption de dépendance des erreurs
library( lmtest)
dwtest(regs)
Durbin-Watson test
DW = 1.0538, p-value = 0.001488
alternative hypothesis : true autocorrelation is greater than 0
On rejette Ho donc il y a autocorélation des erreurs.
Généralement ce problème apparaît sur des données individuelles (données en coupe transversale)
surtout, lorsque les observations représentent des moyennes calculées sur des échantillons de taille
différentes ou bien lorsqu’on répète une même valeur de la variable à expliquer pour des valeurs
différentes d’une variable explicatives, par exemple lors de regroupements par tranches. Il peut aussi
se présenter en cas où les erreurs sont liées aux valeurs privés par une variable explicative, par exemple
la variance de la consommation croit avec le revenu disponible.
a) Test de Goldfeld-Quandt
Le test de Godfeld-Quandt est recommandé lorsqu’on suppose que la variance des erreurs est
fonction d’une variable explicative xk du modèle. Dans ce cas les observations sont ordonnées suivant
25
les valeurs croissantes de xp . Si σ 2 est une fonction croissante de xp alors les observations seront
également ordonnées suivant les valeurs croissantes de la variance des erreurs. Ainsi, ce test n’est
valable que si l’une des variables est la cause de l’hétéroscédasticité et que le nombre d’observations
est important. Ceci autorise à formuler que σi2 = σ 2 x2pi . Le test de Godfeld-Quandt se réalise en trois
étapes :
Étape 1 : on classe les observations de l’échantillon considéré selon l’ordre croissant ou décroissante
de la variable xki .
Étape 2 : On omet de l’échantillon c observations centrales et on divise le reste en deux sous échan-
tillons de même taille n−c
2
.
Étape 3 :On effectue des estimations séparées par MCO sur les deux sous échantillons.
Sous l’hypothèse nulle de l’homoscèdasticité des erreurs le rapport :
SCR2
c ddl2
F = SCR1
∼ F (ddl2 , ddl1 )
ddl1
Yi = a0 + a1 Xi + ui
Pour ce faire, il procède à un test sur 30 véhicules qu ?il regroupe en 6 classes de 5 voitures en
demandant à chaque chef d ?atelier de passer un nombre d ?heures de vérification fixé. Les résultats
sont consignés dans ce tableau
j 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Yj 4 5 6 8 8 6 11 13 15 17 9 13 14 15 21
Xj 4.0 4.0 4.0 4.0 4.0 3.5 3.5 3.5 3.5 3.5 2.0 2.0 2.0 2.0 2.0
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
6 13 16 23 26 11 15 17 22 34 7 21 23 28 38
1.5 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0 1.0 0.5 0.5 0.5 0.5 0.5
Test de Goldfeld-Quandt n’est valable que si l’une des variables est la cause de l’hétéroscédasticité
et que le nombre d’observations est important. Il se compose de trois étapes.
26
Étape 1 : ordonner les observations en fonction des valeurs croissantes ou décroissantes de la variable
à expliquer ou bien de la variable explicative soupçonnée être la source de l’hétéroscédasticité. Cette
opération est déjà effectuée puisque le nombre de défauts constatés est classé (tableau ) en fonction
du nombre d’heures décroissant de vérification.
Étape 2 : omettre C observations centrales. Nous choisissons arbitrairement C observations situées
au centre de l’échantillon. Ces C observations sont exclues de l’analyse. La valeur de C doit être
approximativement égale au quart du nombre d ?observations totales. C = partie entière de (30/4)
=8
Étape 3 : régressions sur les deux sous-échantillons et test.
library(lmtest)
gqtest(regs)
Goldfeld-Quandt test
data : regs
GQ = 6.9463, df1 = 3, df2 = 3, p-value = 0.07286
b) Test de white
Le test de white est basé sur la comparaison des variances estimées des erreurs lorsque le modèle
est estimé par le MCO sous l’hypothèse de l’homoscèdasticité et lorsque le modèle est estimé par le
MCO sous l’hypothèse de l’hétéroscédasticité. Lorsque l’hypothèse H0 est vérifiée, l’écart entre les
deux variances estimé devrait être significativement nul asymptotiqument. En pratique, Le test de
White est fondé sur une relation significative entre le carré du résidu et une ou plusieurs variables
explicatives en niveau et au carré au sein d ?une même équation de régression
27
l ?aide d ?un test de Fisher classique de nullité de coefficients :
H0 : β1 = β2 = ... = βk = βk+1 = 0
nR2 ∼ χ2 (d)
avec d = 2k : le nombre de régression dans l’équation
Si nR2 < χ2α (d), alors l’hypothèse de l’homoscèdasticité est retenue.
Exemple
Soit ici à estimer le modèle :e2j = β0 + β1 xj + β2 x2j + vj
e2j = 136, 02 − 78, 58xj + 11, 98x2j + vj
avec n = 30; R2 = 0, 226; F c = 3, 956
0,05
Test de Fisher F c = 3, 956 > F2;27 = 3, 35
Test LM nR2 = 30 ∗ 0, 22 = 6, 78 > X0,05 2
(2) = 5, 99
Nous sommes, dans les deux cas, amenés à rejeter l’hypothèse H0 pour un seuil de 5 %. Le modèle
est donc hétéroscédastique.
y2=resd 2̂
X2=X 2̂
reg2=lm(y2 ∼ X +X2)
summary(reg2)
où Ω 6= I est la matrice de variance covariance des erreurs qui est symétrique définie positif.
28
2.6.2 L’estimation du modèle linéaire généralisé
L’estimateur de moindres carrés généralisés (MCG) est obtenu en transformant le modèle (??) de
manière à obtenir un modèle de régression linéaire multiple dont les erreurs vérifient les hypothèses
1, 2 et 3 de MCO.
Soit Ω−1 la matrice inverse de Ω symétrique définie positif et qui peut être décomposé telle que
telle que :
Ω−1 = Ω−1/2 Ω−1/2
En multipliant le modèle par Ω−1/2 , on obtient un modèle dont les erreurs vérifie les hypothèses
de MCO.
Soit Y ∗ = Ω−1/2 Y et X ∗ = Ω−1/2 X et ε∗ = Ω−1/2 ε
Le modèle transformé s’écrit :
L’estimateur de MCG est égal à l’estimateur des MCO des paramètres du modèle transformé :
0 0
β̂M CG = σ 2 (X Ω−1 X)−1 X Ω−1 Y
On note que l’application de la méthode de MCG exige de connaître à un constant prés Ω. Lorsque
celle-ci est inconnue, on doit chercher un estimateur de cette matrice de variance covariance Ω̂ dans
ce cas l’estimateur s’appelle moindre carré quasi généralisé (MCQG).
0 0
β̂M CQG = σ 2 (X Ω̂−1 X)−1 X Ω̂−1 Y
Lorsqu’un problème d’autocorrélation d’ordre 1 est détecté, l’estimateur des MCO n’est pas ef-
ficace, il faut appliquer l’estimateur des MCG et par conséquent définir la matrice de variances
covariance des erreurs que s’écrit sous la forme suivant :
1 ρ ρ2 ... ρn−1
1 ρn−2
V ar(ε) = Ω = 1 :
:
1
Afin d’obtenir un estimateur efficace on applique MCG qui consiste à appliquer les MCO sur le
modèle transformé :
Ω−1/2 Y = Ω−1/2 X β + Ω−1/2 ε
29
Avec p
1 − ρ2 0 ... ... 0
−ρ 1 0
2
σ
Ω−1/2
= −ρ 1 :
1 − ρ2
:
1
Par conséquent, pour la première observation, le modèle s’écrit :
p p p p p
1 − ρ2 y1 = β1 1 − ρ2 + β2 1 − ρ2 x11 + ... + βk 1 − ρ2 xk1 + 1 − ρ2 ε1
Lorsque le coefficient d’autocorrélation est inconnu on applique l’estimateur de MCQG. Une méthode
très utilisée est la méthode de Cochrane-Orcutt. Elle procède en plusieurs étapes :
1- Estimer avec la méthode des MCO le modèle initial et évaluer les résidus sur les observations
i = 2, . . . , n.
2- Estimer ρ par ρ̂0 qui est déterminé de manière à minimiser la fonction :
n
X
M in (ei − ρei−1 )2
ρ
i=2
n
P
ei ei−1
On obtient ainsi ρ̂0 = i=2
n . Où à partir de la statistique de Durbin et Waston DW ' 2(1 − ρ̂0 ) :
e2i
P
i=2
3- Appliquer la méthode MCO sur le modèle transformé appelé le modèle quasi-différencie :
yi − ρ̂0 yi−1 = β1 (1 − ρ̂0 ) + β2 (x1i − ρ̂0 x1i−1 ) + ... + βp (xpi − ρ̂0 xpi−1 ) + (1 − ρ̂0 )εi
library(orcutt)
reg2<-cochrane.orcutt(reg)
reg2
Dans le cas de l’héteroscèdasticité des erreurs. La matrice de variance covariance des erreurs peut
s’écrire sous la forme : 2
σ1 0 ... 0
0 0
V ar(ε) = Ω =
: :
0 ... 0 σn2
30
Ce modèle constitue un cas particulier du modèle linéaire généralisé (??) où Ω est une matrice
diagonale (les covariances entre les erreurs sont nulles). La matrice de transformation Ω−1/2 s’écrit :
1/σ1 0 ... 0
0 0
Ω−1/2 =
: :
0 ... 0 1/σn
Plus particulièrement, si le test de l’hétéroscédasticité indique par exemple σi2 = σ 2 x1i on doit
estimer le modèle :
Yi 1 x1i xki εi
√ = β1 √ + β2 √ + ... + βk √ +√
x1i x1i x1i x1i x1i
Enfin lorsque σi2 est inconnue, la méthode de MCQG s’appliquent. Il faut alors dans un premier temps
déterminer un estimateur de la variance des erreurs. Par exemple si : On estime avec la méthode de
MCO l’équation :
e2i = δ1 + δ2 x1i + µi
où µi est une erreur et e2i est le carrée des résidus dans un modèle (1.20) est estimé par le MCO et
ê2i = δ̂1 + δ̂2 x1i est un estimateur de e2i . Finalement, on calcule l’estimateur MCQG en estimant par
les MCO :
Yi 1 x1i xpi εi
= β1 + β2 + ... + βp +
êi êi êi êi êi
31
Chapitre 3
Jusqu’à maintenant, nous avons spécifié des modèles où toutes les variables sont exprimées à la
même période. Cependant, la théorie économique postule couramment non pas des effets synchrones
mais des effets retardés.
Ces modèles à décalages peuvent inclure comme variables retardées aussi bien la variable endogène
que les variables exogènes. Nous étudierons, tout d’abord en section1., les modèles autorégressifs
dont la variable retardée est la variable endogène. Puis, en section 2., une autre classe de modèle
est envisagée : les modèles à retards échelonnés. Dans cette spécification, les variables exogènes
apparaissent à plusieurs décalages. Enfin, en section 3, nous présentons deux exemples de modèles
dynamiques.
Dans ce modèle, l’hypothèse d’indépendance entre les variables explicatives et l’erreur, n’est plus sa-
tisfaite car les variables yt−1 , yt−2 , ..., yt−h qui dépendent de t−1 , t−2 , ..., t−h sont aléatoires puisque
yt+1 est fonction de yt qui dépend de t , soit : E(yt+1 t ) 6= 0
Or nous devons considérer les processus explosifs que ce soit d’une façon monotone ou oscillatoire ?
comme des phénomènes rares en économie (les variables sont limitées généralement dans leur crois-
sance), c’est pourquoi il convient d’ajouter l’hypothèse de stabilité du processus.
32
Cas particulier : le modèle autorégressif d ?ordre 1
Le modèle général spécifié ci-dessus est rarement utilisé, le plus souvent nous nous limitons à des
processus autorégressifs d ?ordre 1, de la forme :
yt = b0 + b1 yt−1 + a1 x1t + t
Ce modèle est stable si |b1 | < 1 et explosif si |b1 | > 1 .
a) Test d’autocorrélation des erreurs Dans le cas d’un modèle autorégressif d ?ordre 1, le test
classique de Durbin et Watson sous-estime le risque d ?autocorrélation, un nouveau test d’autocor-
rélation des erreurs doit alors être utilisé. La statistique utilisée est la suivante :
Il faut utiliser la statistique h de Durbin pour tester l’auto corrélation d’ordre 1 des erreurs :
DW n
r
h = (1 − )
2 1 − nV (b̂1 )
L’hypothèse nulle ρ = 0 est acceptée si |h| ≤ 1.645 pour un risque de première espèce de 5%. Le
test est unilatéral, l’hypothèse alternative est ρ > 0 ou ρ < 0
Si nV (b̂1 ) prend une valeur supérieur à 1, on ne peut pas calculer la statistique h de Durbin. On
peut procéder de la manière suivante :
1-Estimer le modèle dynamique avec la méthode des MCO.
2-Estimer le modèle suivant :
ei = α1 + α2 x2,i + α3 ei−1 + εi i = 1, .., n
où ei est le résidu de MCO du modèle dynamique.
3-Si le coefficient de la variable ei−1 est différent de 0 on peut rejeter l’hypothèse nulle de non auto
corrélation.
b) Estimation si les erreurs ne sont pas corrélées Dans le cas d’absence d’ auto corrélation
des erreurs, les estimateurs des MCO appliqués au modèle convergent asymptotiquement vers les
valeurs vraies des paramètres et ont une variance minimale parmi tous les estimateurs convergents.
Cependant, pour les petits échantillons, lors de l’estimation d’un modèle autorégressif d ?ordre h
, les résultats asymptotiques sont alors très approximatifs, car le nombre de périodes d’estimation
est de n - h. De plus, les problèmes de colinéarité entre les variables explicatives décalées interdisent
pratiquement d’utiliser les MCO.
En résumé, pour l’estimation des modèles autorégressifs à erreurs indépendantes, l’application de
la méthode classique des MCO est licite si le nombre d’observations est suffisant (souvent dans la
pratique n > 15).
33
c) Estimation en cas d’autocorrélation des erreurs Le modèle s’écrit alors :
Avec t = ρt−1 + vt
Il s’agit d ?un modèle à autocorrélation des erreurs, nous avons montré lors du chapitre précédant
que la transformation en quasi-différences premières peut lever l’autocorrélation des erreurs d ?ordre
1. Nous pouvons citer différentes méthodes d’estimation.
a) Régression sur les différences premières :On procède à deux régressions par la méthode des MCO :
la première sur le modèle spécifié avec des variables non transformées et la seconde sur le modèle
spécifié avec des variables transformées en différences premières. Si les deux résultats d’estimation
sont proches, nous pouvons conclure que les paramètres du modèle sont suffisamment bien déterminés
sans rechercher la liaison éventuelle des erreurs.
b) Correction de l’autocorrélation des erreurs par régression sur les quasidifférences premières, nous
procédons à la correction de l’autocorrélation des erreurs : estimation directe de ρ , Hildreth-Lu,
Cochrane-Orcutt.
Application
Un économétre désire tester la relation entre les prix officiels (PO) de la tonne de café et les prix
réellement appliqués à l’exportation (PE) par les pays producteurs. Il se propose d’estimer la relation :
P Ot = a0 + a1 P Ot−1 + a2 P Et + t
dans laquelle le prix officiel est fonction de manière instantanée du prix observé et s’ajuste avec un
retard d’un an au prix officiel. Il dispose des données, sur 16 ans, présentées dans ce tableau :
L’équation est estimé sur 15 périodes (n -1) car dans le modèle figure une variable retardée d ?une
période. Les résultats de l’estimation sont les suivants :
34
La statistique de Durbin et Watson laisse augurer d’une autocorrélation des erreurs, ce qui est
confirmé par le « h » de Durbin :
r
DW n 15
r
h = (1 − ) = 0.687 ∗ = 2.68 > 1.96
2 1 − nV (b̂1 ) 1 − 15 ∗ 0.0013
En général, l’effet de la variable explicative s’estompe avec le temps a0 > a1 > a2 > ... > ah
Ce modèle peut se simplifier dans son écriture en utilisant un opérateur décalage définit par : Dxt =
xt−1 et en général : Di xt = xt−i , soit :
h h
" h
#
X X X
yt = aj xt−j + b0 + εt = aj D j x t + b 0 + ε t = aj D j x t + b 0 + ε t
j=0 j=0 j=0
yt = A (D) xt + b0 + εt
Où A (D) est polynôme de degré h, tel que :
A (D) = a0 + a1 D1 + a2 D2 + ... + ah Dh
Le nombre de retards, h, peut être fini ou infini. Cependant, la somme des coefficients aj tend vers
une limite finie, si non yi serait un processus explosif.
Le polynômeA (1) = a0 +a1 +a2 +...+ah permet de mesurer, à long terme, l’impact de la modification
de la variable explicative xt d’une quantité ∆x sur la valeur deyt . En effet, les coefficients aj représente
les multiplicateurs instantanés et leurs sommes le multiplicateur cumulé.
35
3.2.2 Retard moyen
La notion de retard moyen permet d’appréhender le laps de temps qui correspond à la valeur
moyenne des coefficients aj . Il est égale à la moyenne pondérée des coefficients, soit :
Ph
j=0 jaj a1 + 2a2 + ... + hah A0 (1)
D̄ = Ph = =
j=0 aj a0 + a1 + a2 + ... + ah A (1)
L’estimation des paramètres du modèle (1) soulève deux types de difficultés :
1. En pratique, la valeur h qui mesure l’importance du nombre de retards est rarement connue,
même si nous en connaissant l’ordre de grandeur.
2. Le second problème résulte de la colinéarité entre les variables exogènes décalées. En effet,
lorsque le nombre de retards est important la colinéarité entre les variables explicatives décalées
risque d’entraîner une imprécision dans l’estimation des coefficients. C’est pourquoi, le plus
souvent, nous émettons des hypothèses sur la forme de retards.
Ces deux points sont abordés ci-dessous :
a)Test de Fisher
Le test le plus naturel est celui de Fisher dans lequel nous testons l’hypothèse de la nullité des
coefficients de régression pour les retards supérieurs à h.
La formulation des hypothèses est la suivante, lorsque l’on teste, d’une manière descendante, une
valeur de h comprise entre 0 et M 0 < h < M .
H01 : h = M − 1 → aM = 0 H11 : h = M → aM 6= 0
H02 : h = M − 2 → aM −1 = 0 2
H1 : h = M − 1 → aM −1 6= 0
.. ..
. .
H0i : h = M − i → aM −i+1 = 0 H1i : h = M − i + 1 → aM −i+1 6= 0
36
b) Critère de Akaike (AIC)
Une autre méthode consiste à retenir comme valeur de h celle qui minimise la fonction de Akaike qui
est donnée par :
SCRh 2h
AIC (h) = Ln +
n n
avec SCRh = Sommes des Carrés des Résidus pour le modèle à h retards
n = Nombre d’observations
Ln = Logarithme népérien
Enfin, une méthode très proche de la précédente consiste à retenir la valeur de h qui minimise la
fonction de Schwarz :
SCRh hLnn
SC (h) = Ln +
n n
Exemple : La théorie économique postule que les dépenses d’investissement (notées yt ) peuvent
être expliquées par les profits passés (notés xt ). Le modèle prend la forme d’un modèle à retards
échelonnés.
N’étant pas certains de la spécification exacte du modèle, nous désirons tout d’abord rechercher le
nombre de décalages trimestriels qui semblent avoir un effet sur les dépenses d’investissement. Puis
après avoir déterminé le nombre de retards, nous calculerons le délai moyen.
Nous procédons aux estimations des modèles sur 34 trimestres, puisque le décalage maximum testé
est de 10, on élimine les 10 premières observations. Nous calculons les trois critères : Akaike, Schwarz
et Fisher à l’aide d’un programme en langage R.
37
tab2 = read.csv("retard echolonne.csv", header=TRUE, sep=" ;", dec=",")
attach(tab2)
y=cbind(Y,X)
y.ts =ts(y, start=c(1990,1), end=c(2000,4),frequency=4)
Y=y.ts[,1]
X=y.ts[,2]
library(dynlm)
aics =rep(0,10)
bics = rep(0,10)
for (i in 1 :10)
ari = dynlm(Y∼L(X,1 :i), start=i)
aics[i] = AIC(ari)
bics[i] = BIC(ari)
library(knitr) # for kable()
tbl = data.frame(rbind(aics, bics))
names(tbl) = c("1","2","3","4","5")
row.names(tbl) =c("AIC","BIC")
kable(tbl, digits=1, align=’c’,caption="Lag order selection for an AR model")
∞ ∞
" ∞
#
X X X
yt = aj xt−j + bo + εt = aj D j x t + b o + ε = aj D j x t + b 0 + ε t
j=0 j=o j=o
Afin de ce ramener à un nombre fini de paramètres a estimé nous devons postuler une forme
particulière que peut prendre la succession des coefficients ao , a1 , a2 ...
Nous pouvons admettre plusieurs types de spécifications, de modèles particuliers –les plus utilisés-
font l’objet d’une présentation :
1. Une des croissances géométriques des effets de la variable exogènes avec le temps ;
2. Une croissance suivie d’une décroissance.
a1 = λa0
a2 = λ2 a0
38
Et en général ai = λi ao avec o < λ < 1
Soit
yt = b0 + a0 xt λDxt + λ2 D2 xt + .. + λi Di xt + εt
y t = b 0 + a0 x t S + ε t
Or S = (1 + λD1 + λ2 D2 + ...) = (1 − λD)−1 (somme d’une progression géométrique), nous avons
donc
(1 − λD) yt = a0 xt + (1 − λ) b0 + (1 − λD) εt
yt = λyt−1 + a0 xt + (1 − λ) b0 + εt − λεt−1
Qui est un modèle autorégressif à erreurs auto corrélées dont la procédure d’estimation a été exposée
en I.
Cette transformation qui consiste à passer d’un modèle à retards échelonnés à un modèle autorégressif-
est habituellement appelée « transformation de Koyck ». Nous pouvons noter que nous passons d’un
modèle à retards échelonnés, difficile à estimer par l’abondance des paramètres, mais à erreurs ré-
pondant aux hypothèses classiques, à un modèle autorégressif, simple dans sa spécification, mais à
erreurs auto corrélées.
yt = A (D) xt + b0 + εt
39
Une démonstration analogue à la précédente permet d’écrire le modèle [14]
B (D) yt = B (D) A (D) xt + B (D) b0 + B (D) εt avec B (D) = A (D)−1
avec
- Pour r = 0 → B (D) = (1 − λD) /a0 , soit le modèle précédent
- Pour r = 1 → B (D) = (1 − λD)2 /a0
Soit
yt = 2λyt−1 − λ2 yt−2 + a0 xt + 1 − 2λ + λ2 b0 + vt
Avec vt → AR (2)
- Pour r = 2 → B (D) = (1 − λD)3 /a0
Soit :
Avec v1 → AR (3)
Afin de déterminer les valeurs du paramètre r, Maddala et Rao (1971) suggèrent l’utilisation d’une
procédure de balayage dont la fonction objective à maximiser est le coefficient de détermination
corrigé R̄2 .
– Afin de faire face à une augmentation de la demande, une entreprise cherche à se doter de
moyens de production supplémentaires. Cependant cet investissement ne peut pas être réalisé
immédiatement et demande un certain temps d’ajustement.
– Lors de l’augmentation des prix du pétrole, les consommateurs et les institutions réagissent
avec retard à cette modification des prix. Il en résulte qu ?à court terme le comportement des
consommateurs ne reflète pas l’ajustement qui s’est produit à long terme
La formulation de ce modèle est la suivante : il est nécessaire de distinguer entre la valeur désirée,
de la variable à expliquer ytD , et la valeur vraie de cette même variable (yt ) . Le niveau désiré de la
variable à expliquer est fonction de la variable explicative xt tel que :
ytD = a0 + a1 xt + t
40
Une modification de la valeur de xt entraîne une modification de la valeur désirée de la variable
à expliquer. Cependant, malgré une spécification simple, nous ne pouvons pas appliquer les MCO
puisque nous n’avons aucune mesure de ytD . En revanche, nous connaissons les valeurs de yt et nous
pouvons spécifier une relation entre ytD et yt :
yt = a0 + a1 xpt + εt
Où xpt est la valeur prévue de la variable explicative xt .
Par exemple la production d’une entreprise est fonction des valeurs prévues des ventes. Comme pour
le modèle précédent, nous sommes confrontés à l’absence de mesure de la variable xpt . Pour lever
cette difficulté, nous devons poser une hypothèse concernant la formulation de la variable xpt : celle
des anticipations adaptatives.
Cette relation s’écrit :
41
i=∞
X
y t = a0 + a1 λ (1 − λ)i xt−1 + εt
i=0
Qui est un modèle à retards échelonnées dans une transformation de Koyck permet de mettre sous
une forme autorégressive :
42
Chapitre 4
La Cointégration
43
4.2 Test de cointégration : Test de Engel et Granger
Ce test est construit en deux étapes :
Etape 1 : tester l’ordre d’intégration des variables
Une condition nécessaire de cointégration est que les séries doivent être intégrées de même ordre. Si
les séries ne sont pas intégrées de même ordre, elles ne peuvent pas être cointégrées.
Étape 2 : estimation de la relation de long terme
Si la condition nécessaire est vérifiée, on estime par les MCO la relation de long terme entre les
variables :yt = a1 xt + a0 + εt
Pour que la relation de cointégration soit acceptée, le résidu et issu de cette régression doit être
stationnaire : et = yt a1 xt − a0 .
La stationnarité du résidu est testée à l’aide des tests Dickey-Fuller et Dickey-Fuller Augmenté.
En effet, le test porte sur les résidus estimés à partir de la relation statique et non pas sur les « vrais
» résidus de la relation de cointégration. Mackinnon (1991) a donc simulé des tables qui dépendent
du nombre d’observations et du nombre de variables explicatives figurant dans la relation statistique.
Si le résidu est stationnaire nous pouvons alors estimer le modèle à correction d’erreur.
yt = α0 + α1 yt−1 + α2 xt + α3 xt−1 + vt
La dynamique de long terme s’exprime de la manière suivante :
yt = axt + b + εt
Car à long terme, on a yt−1 = yt ,xt−1 = xt et la dynamique de court terme devient à long terme :
yt = α0 + α1 yt + α2 xt + α3 xt + vt
D’où (1 − α1 )y t = (α2 + α3 )xt + α0 + vt
yt = axt + b + εt
Avec = α1−α
2 +α3
1
α0
, b = 1−α 1
vt
,εt = 1−α 1
.
Le MCE s’obtient à partir de la dynamique de court terme :
yt = α0 + α1 yt−1 + α2 xt + α3 xt−1 + vt
44
4y t = (α1 − 1) yt−1 + α2 (xt − xt−1 ) + α0 + (α2 + α3 )xt−1 + vt
α2 + α3 α0
4y t = − (1 − α1 ) yt−1 − xt−1 − + α2 4xt + vt
1 − α1 1 − α1
yt , xt ∼> CI(1, 1)
Nous pouvons estimer le MCE.
- Étape 1 :estimation par les MCO de la relation de long terme :
yt = axt + b + εt
– Étape 2 : estimation par les MCO de la relation du modèle dynamique (court terme) :
45
Chapitre 5
Parfois la spécification nécessite l’écriture d’équations multiples reliées entre elles au travers des
variables figurants dans plusieurs équations.Nous ne pouvons pas, sauf cas particuliers, utiliser la
méthode des MCO équation par équation comme si chacune de ces équation était indépendante les
unes des autres. Nous envisageons la présentation matricielle d’un modèle à équations simultanées
et à son écriture sous forme réduite (les variables endogènes sont exprimées en fonction des seules
variables exogènes).Le problème de l’identification, c’est-à-dire la possibilité d’estimation du mo-
dèle sous forme réduite. Enfin, nous présentons les méthodes d’estimation spécifique des modèles à
équations simultanées.
46
Ce système d’équations multiples, spécifiés par l’économiste, qui traduit directement les relations
entre les variables, s’appelle : le systèmes d’équations structurelles. Ce modèle comporte trois équa-
tions dont une identité [E3]. En effet, dans la relation [E3], il n’y a aucun coefficient à estimer et par
conséquent pas de termes aléatoires. L’équation [E1] est une fonction de consommation et l’équation
[E2] est relative à l’investissement.
Ce système contient trois variables endogènes Ct ,It , Yt est une variable exogène Yt−1 . Nous re-
marquons, par exemple, que la variable Yt apparaît comme variable explicative en [E1], ce qui est
contraire à son statut de variable endogène. Pour lever ce problème, nous allons exprimé les trois
variables endogènes (Ct ,It ,Yt ) en fonction de la seule variable exogène (Yt−1 ).
En substituant [E3] dans [E1], nous obtenons :
a0 + a1 b 0 a1 b 1 a1 ε2t + ε1t
Ct = + Yt−1 +
1 − a1 1 − a1 1 − a1
Il en résulte que :
h i
a0 +a1 b0 a1 b1 a1 ε2t +ε1t
Yt = Ct + It = 1−a1
+ b0 + 1−a1
+ b1 Yt−1 + 1−a1
+ ε2t
a0 +b0 b1 ε2t +ε1t
Yt = 1−a1
+ Y
1−a1 t−1
+ 1−a1
Les équations structurelles sont alors équivalente aux équations réduites (les variables endogènes sont
exprimés en fonction des seules variables exogènes) :
Ct = a01−a
+a1 b0
1
a1 b1
+ 1−a 1
Yt−1 + a1 ε1−a
2t +ε1t
1
[E4]
a0 +b0 b1 ε2t +ε1t
Yt = 1−a1 + 1−a1 Yt−1 + 1−a1 [E5]
It = b0 + b1 Yt−1 + ε2t [E6]
L’équation [E5] indique que la variable Yt est fonction de ε1t et par conséquent E (Yt ε1t ) 6= 0. Il
en résulte que, dans l’équation [E1], l’hypothèse d’indépendance entre « variables explicatives Yt et
l’erreur ε1t n’est pas respecté et l’application des MCO sur le modèle [E1] conduit à des estimateurs
biaisés et non convergents.
En revanche, l’utilisation des MCO sur les équations réduites est licite puisque la variable Yt−1 est
indépendante de ε1t et ε2t .
Il est noté que la forme réduite permet de mesurer l’effet total, direct et indirect, d’une modification
de la variable exogène Yt−1 sur les variables endogènes. Par exemple, le schéma suivant illustre la
cascade des causalités :
[E2] [E3] [E1]
47
Yt−1 → It → Yt → Ct ce qui est résumé dans [E4]
b11 y1t + b12 y2t + ... + b1g ygt + c11 x1t + c12 x2t + ... + c1k xkt = ε1t
b21 y1t + b22 y2t + ... + b2g ygt + c21 x1t + c22 x2t + ... + c2k xkt = ε2t
...
bg1 y1t + bg2 y2t + ... + bgg ygt + cg1 x1t + cg2 x2t + ... + cgk xkt = εgt
Soit, sous forme matricielle :
B Y + C X = ε
(g,g) (g,1) (g,k) (k,1) (g,1)
Bien entendu, dans chaque équation, quelques coefficients sont nuls et la variable dont le coefficient
est égale à 1 et la variable dépendante. Si une équation à tous ses coefficients, non nuls, égaux à 1 et
ne comportent pas de terme aléatoire, cela signifie qu’il s’agit d’une identité (aucun coefficient n’est
à estimer).
Si la matrice B est régulière, nous passons de la forme structurelle à la forme réduite en exprimant
Y en fonction de X, soit :
Y = −B −1 CX + B −1 ε
Et nous pouvons appliquer les MCO. En effet, les erreurs B −1 ε sont indépendantes de x. Si, sur
le plan de la représentation la formalisation est simple, sa mise en oeuvre pratique est plus complexe.
En effet, la connaissance de g × k élément de la matrice (−B −1 C)ne permet pas de déterminer –
d’identifier- la matrice B (composé de g × g éléments) ainsi que la matrice C (composé de g × k
éléments). Nous sommes alors en présence de g × k équations à (g × g) + (g × k) inconnues qui,
sont restrictions supplémentaires, s’avèrent impossibles à résoudre. A titre d’illustration le modèle
introductif peut s’écrire sous forme matricielle :
Ct 1 0 −a1 " # −a0 0
U
Y = It ; B = 0 1 0 ;X = ; C = −b0 −b1
Yt−1
Yt −1 −1 1 0 0
ε1t
ε = ε2t (U est le vecteur unité). Nous allons donc présenter dans la section suivante des règles
0
simples permettant, à partir de restriction sur les paramètres, de déterminer des conditions d’iden-
tification.
48
5.2 Le problème de l’identification
5.2.1 Restrictions sur les coefficients
Il y a une restriction sur un coefficient de la forme structurelle, chaque fois qu’un paramètre est
contraint –par l’écriture du modèle- à être égale à une valeur déterminée. Nous distinguons deux
types de restrictions.
a)Restrictions d’exclusion
Nous pouvons considérer que chaque fois qu’une variable endogène ou exogène n’apparaît pas dans
une équation structurelle, cela revient à l’affecter d’un coefficient nul. Par exemple, dans notre modèle
introductif, la variable It ne figure pas dans l’équation [E1], son coefficient est donc nul : dans la
matrice B, l’élément de la première ligne est de la deuxième colonne égale à 0.
b)Restrictions linéaires
Certaines spécifications de modèles imposent que des variables soient affectées d’un coefficient iden-
tique, il s’agit là encore de restriction à priori sur les paramètres du modèle.
49
g − 1 > g − g 0 + k − k 0 →l’équation est sous- identifiée
g − 1 = g − g 0 + k − k 0 → l’équation est juste identifiée
g − 1 < g − g 0 + k − k 0 → l’équation est sur- identifiée
Ce qui peut se résumer ainsi : pour qu’une équation ne soit pas sous identifiée, le nombre de variables
exclues de l’équation doit être au moins égale au nombre d’équations du modèle moins 1.
Lorsque nous avons r restrictions, autres que celles d’exclusion, concernant les paramètres d’une
équation (égalité des deux coefficients, par exemple), les conditions précédentes deviennent :
g − 1 > g − g 0 + k − k + r → l’équation est sous identifiée
g − 1 = g − g 0 + k − k 0 + r → l’équation est juste identifiée
g − 1 < g − g 0 + k − k 0 + r → l’équation est sur- identifiée
Ces conditions –nécessaires- sont appelées conditions d’ordre d’identifiabilité.
Il convient de vérifier des conditions suffisantes, qualifiées de conditions de rang qui dans la pratique,
se révèle difficile, voir parfois impossible de mettre en oeuvre.
50
plus simple à mettre en oeuvre et qui fournit les mêmes résultats que les MCI pour les équations
justes identifiées.
b11 y1t + b12 y2t + ... + b1g ygt + c11 x1t + c12 x2t + ... + c1k xkt = ε1t
b21 y1t + b22 y2t + ... + b2g ygt + c21 x1t + c22 x2t + ... + c2k xkt = ε2t
...
bg1 y1t + bg2 y2t + ... + bgg ygt + cg1 x1t + cg2 x2t + ... + cgk xkt = εgt
La première étape consiste à effectuer une régression de chacune des variables endogènes sur les
variables exogènes :
y1t = β12 ŷ2t + ... + β1g ŷgt + c11 x1t + c12 x2t + ... + c1k xkt + ε1t
y2t = β21 ŷ1t + ... + β2g ŷgt + c21 x1t + c22 x2t + ... + c2k xkt + ε2t
...
ygt = βg1 ŷ1t + ... + βgg ŷgt + cg1 x1t + cg2 x2t + ... + cgk xkt + εgt
Les propriétés de l’estimateur des DMC sont identiques, de manière asymptotique, à celle d’un esti-
mateur classique ; c’est-à-dire que pour les petits échantillons les estimations des paramètres peuvent
être biaisées. Il a noté que l’estimateur des DMC peut s’interpréter comme étant un estimateur des
variables instrumentales, les variables exogènes des autres équations étant les instruments.
51
Profit : profits du secteur privé
Capital : stock de capital de l’année précédente
Govtwage : salaires versés par le gouvernement
Govtexp : dépenses de l’Etat en dehors des salaires
Taxes : taxes
Residuals :
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.8930 -0.6163 -0.2459 0.0000 0.8854 2.0030
52
Instruments : ∼G + T + Wg + I(Year - 1931) + K.lag + P.lag + X.lag
Residuals :
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.2910 -0.8069 0.1423 0.0000 0.8601 1.7960
Residuals :
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.2940 -0.4734 0.0145 0.0000 0.4487 1.1960
53
Chapitre 6
L’analyse économétrique se base sur la spécification d’une relation théorique. On commence l’es-
timation du modèle qui est supposé correctement spécifié, avec des hypothèses optimistes à partir de
lesquelles on peut obtenir des mesures précises de tous les paramètres des variables. Si les conditions
idéales sont satisfaites à chaque étape, l’analyse est probablement routinière. Malheureusement, ces
conditions ne sont que rarement satisfaisantes. Ainsi, des difficultés que l’on peut rencontrer sont :
- Les propriétés des termes aléatoires supposées stochastiques peuvent être violées.
- Des caractéristiques spécifiques sont corrélées avec des variables explicatives.
- Des effets individuels existent dans le modèle.
Vue leurs doubles dimensions, les données de panel peuvent s’interpréter comme solution de ces
problèmes : c’est le domaine le plus actif et le plus navetteur de l’économétrie, car les données de
panel constituent un cadre d’analyse favorable au développement des techniques d’estimation et des
résultats théoriques.
La disponibilité d’échantillon de données de panel est considérablement développée d’une part
du fait de la prise de conscience de l’intérêt de ce type de données mais aussi bien grâce au progrès
informatique. Dans la pratique les données de panel permettent d’étudier des questions impossibles
à traiter en coupe transversale ou en série temporelle. Le développement de la théorie économique
vise à une meilleure représentation de la dynamique des comportements des agents et à la prise en
compte de leur hétérogénéité. Les données de panel, part leur double dimension se révèlent donc
particulièrement adoptées dés lors que l’on souhaite estimer ces modèles et tester ces théories. Il est
par ailleurs remarquable que pendant très longtemps, les développements de l’économétrie de données
de panel c’est surtout centré sur les problèmes posés par l’estimation de modèles microéconomiques.
Ce chapitre est consacré à présenter les méthodes d’estimations des modèles ajustées aux données
de panel. Nous illustrerons les méthodes proposées par une application qui consiste à examiner
l’efficacité des services dans le transport aérien. En effet, nous analyserons la fonction du coût de
production de quelques firmes de transport américain.
Dans cette étude, on essayera d’illustrer les modèles des données de panel en utilisant des tests
appropriés qui permettent de distinguer le modèle à effets fixes et les modèle à effets aléatoires.
L’objectif de notre étude empirique est d’examiner l’efficacité des services du transport aérien en
analysant le facteur du coût de production d’un panel constitué de 6 firmes de transport aérien
54
Américain pour la période 1974-1985. En réalité, le non disponibilité des données de panel et surtout
en ce qui concernée de firmes tunisienne aussi bien que le tems nous a exigé le choix de ces firmes et
de ces périodes. Les paramètres caractérisent la fonction du coût de production sont en fonction des
variables exogène suivantes.
- La production est mesurée en terme de revenue (passage par nul).
- Le taux de remplissage moyen des sièges des avions représente la capacité d’utilisation.
- Le prix de carburant.
L’originalité de cette étude tient au faite qu’on cherche à spécifier, l’effet de ces variables exogène
comme support de l’efficacité de service dans le transport aérien.
N X
T N N X
T
(Ȳi − Ȳ¯ )2 +
X X X
(Yit − Ȳ )2 = T (Yit − Ȳi )2
i=1 1=t i=1 i=1 t=1
T N
P ¯ 2
i=1 (Ȳi − Ȳ ) : traduit les différences permanentes entre les individus.
PN PT 2
i=1 t=1 (Yit − Ȳi ) : traduit des écarts entre la situation des individus a chaque date.
55
présentent un double avantage d’une part ils conduisent dans les procédures de confrontation des
théories (microéconomiques). Au fait observées aune meilleure adéquation entre le niveau d’analyse
de modèle théorique et celui des observations statistiques, d’autre part il permet d’éviter certaines
difficultés liées a l’agrégation.
Ce modèle décrit un système de N équations (une équation pour chaque individu i). Chaque
équation étant composé de T observations. Si les aléas sont non auto corrélés et si les aléas inter
équation ont une covariance nulle, ce modèle peut être estimer par la méthode de MCO. Si les aléas
sont auto corrélées, il faut appliquer les MCG. Si les aléas inter équation ont un covariance non nulle
56
il faut appliquer de méthode de SUR.
Si dans le modèle précédant ai = aet bi = b et si var(εit ) = σ 2 etcov(εit , εjt ) = 0alors le modèle est
un modèle à une équation et il doit être estimer par la méthode MCO . Si non l’équation doit être
estimée par la méthode MCG
Si ai 6= a et bi = b, ∀i = 1, . . . , N alors le modèle est un modèle à une équation avec des effets
individuels. De manière générale le modèle a effet individuel s’écrit :
b0it = b0 + ai + dt
Par exemple, On suppose que Yi représente la production d’un pays i a la date t et Xit représente
les divers variables explicatives (travail, capital, le niveau d’éducation). D’après ce modèle, à une date
donnée, deux pays ayant les mêmes caractéristiques observées doivent avoir, à une constant près, le
même niveau de production.
(En espérance :E(Yit /X1it , . . . , Xkit ) = b0 + ai + dt + kk=1 bk Xkit )
P
Une différence de produit entre ces deux pays si elle existe, est alors liées a des spécificités inobser-
vables stables dans le temps, dont l’effet est mesuré par le coefficientai . Par exemple cette différence
est du aux facteurs socioculturels ce différence est difficile a quantifiés, mais dont l’effet nous intéressé
néanmoins.
Dans ce modèle on suppose que l’hétérogénéité des comportements est située uniquement dans la
dimension individuelle.
XK
Yit = b0 + ai + bk Xkit + εit
k=1
57
ai : L’effet fixe individuel
Pour rendre plus explicite le fait que ces effets fixes individuels sont de coefficients à estimer on peut
introduire dans l’écriture du modèle les variables indicatrices Zj associées a chacun de ce coefficients
aj :
X K
X
Yit = b0 Z0it + aj Zjit + bk Xkit + εit
k=1
(
1 ∀t, si j = i
Z0it = 1, ∀(i, t) et Zjit =
0∀t, si j 6= i
Le théorème de Frisch -Waugh montre comment l’estimation de coefficient d’un modèle économé-
trique peut se décomposer en celle, séquentielle, de deux sous ensemble de ses coefficients. Supposant
ainsi que l’on souhaite estimer par les MCO le modèle sans constante,
Mz Y = Mz b + Mz ε
58
1 0
. . ···
ε11
.
Z(N T,N ) = 0
1 0 .ε(N T,1) =
ε1T
..
.
.
0
. 1
εN t
0 . 1
En effet la matrice Z peut être réécrite de façon compact comme Z = IN ⊗ eT , ou IN est la
matrice identité d’ordre N et eT est un vecteur colonne unitaire de dimension T et ou⊗ représente
le produit de Kronecker des matrices. Par conséquent :
0
Z 0 Z = (IN ⊗ eT )(IN ⊗ et ) = IN ⊗ T
et (Z 0 Z)−1 = IN ⊗ T1
Ainsi
b̂intra−in = (X 0 WN X)−1 X 0 WN y
V ar(b̂intra−in) ) = σ̂w2 (X 0 WN X)−1
Sous l’hypothèse que les perturbations ε suivent une loi normale Avec b̂inta−in ∼ N (b, σ̂w (X 0 WN X)−1
Modèle simple :
59
Modèle à effets fixes individuels :
fixef(air.fe)
1 2 3 4 5 6
9.705942 9.664706 9.497021 9.890498 9.729997 & 9.793004
Le Tableau présente les résultats de l’estimation du modèle a effets fixes par les estimateurs des
MCO pour l’échantillon de 6 firmes pendant les 15 périodes. A partir de ce tableau, on peut tirer les
constatations suivantes :
-La qualité de l’ajustement R2 indique que 99% des variations de la variable endogène sont expliques
par le modèle. Par ailleurs le coefficient associé à la constante est éliminé puisque le modèle est à
60
effet individuel (effet de firme).
-Conformément à la théorie, la production agit positivement sur le coût. Un accroissement de 1%de
la production induit a un accroissement de 91% du coût. En fin le taux du remplissage à un impact
négatif sur le coût, une augmentation du taux de remplissage induit une diminution du coût.
-L’estimateur des MCO des paramètres du modèle à effet fixe individuel est identique à l’estimateur
de Within qui applique l’estimateur de MCO dont les données sont centrées par rapport a la moyenne
individuel.
Nous avons supposé l’existence d’une hétérogénéité inter – individuelle des comportements, hé-
térogénéité prise en compte par l’adjonction d’effets individuels fixes aux « véritables » variables
explicatives. Il parait naturel de tester l’existence d’une telle hétérogénéité. Ceci revient à discrimi-
ner entre le modèle :
XK
yit = b0 + bk xkit + εit (1)
k=1
et le modèle avec effets individuels fixes :
K
X
yit = a∗i + bk xkit + εit (2)
k=1
61
6.3.2 Le modèle à effets aléatoires individuels
La perturbation d’un modèle économétrique représente l’ensemble des influences des variables non
prise en compte dans le modèle. Ces variables peuvent être absents par décision de l’économètre (qui
considère, éventuellement à tort, qu’elles n’ont pas d’impact sur la variable à expliquée) ou malgré
lui s’il ne dispose pas des observations statistiques sur ces variables. Par exemple la production
des entreprises dépend des quantités de facteurs utilisés (capital, travail) mai aussi de qualité de ces
facteurs (technologie incorporée dans les équipements aux niveaux de qualification de main d’oeuvre)
a) Spécification du modèle
k
X
yit = b0 + bk xit + εit ∀i = 1, .., N et t = 1, .., T
k=1
La perturbation εit est donc composée de deux éléments, d’où la dénomination du modèle repré-
sente l’effet individuel ;ui rendant compte de l’influence sur y des variables non prises en compte wit
représente l’influence des autres variables omises variant d’un individu à l’autre. Ainsi ui e twit sont
des variables aléatoires d’espérance nulle non auto crollées mais pouvant être corrélées elle, et en plus
les variables explicatives du modèle sont supposées strictement exogènes et les deux composantes de
la perturbation vérifient l’hypothèse suivant
cov(yit ; yit0 /x1it , x2it , dots, xkit ) = cov(εit , εi0 t0 /x1it , x2it , dots, xkit )
2 2 0 0
σu + σw si i = i et t = t
= σu2 si i = i et t 6= t0
0
σu2 + σw2 · · · σu2
σu2 σu2
...
= σw2 IT +
Soit A la matrice de variance covariance des aléas A(T,T ) = .. ..
. .
σu2 σu2 σu2 + σw2
σu2 JT Ou IT est la matrice identité d’ordre T et JT est une matrice unitaire (T, T )
62
b) Estimation du modèle à effets aléatoires individuels
En empilant l’ensemble des observations individu par individu, on peut écrire matriciellement le
modèle :
b̂ = (X 0 WN X + θX 0 BN X)−1 (X 0 WN y + θX 0 BN y)
b̂(λ) = X 0 WN X + λX 0 BN X)−1 (X 0 WN Y + λX 0 BN Y )
Est égale
2
σW
λ= 2
=θ
σW + T σu2
Par ailleurs, la mise en ouvre de cet estimateur est assez simple. En effet comme tout estimateur
de MCG il peut s’obtenir l’estimateur de MCO appliqué au modèle transformé suivant :
−1 −1 −1
A( 2
)
Y = A( 2
)
Xb + A 2 ε
L’estimateur de MCG ne peut être calculé, il faut dans un premier temps estimer les variancesσu2 etσw2
pour en déduire une estimation deθ. La substitution de cette estimation θ̂ a la vraie valeur (inconnue)
deθ permet de calculer ce qu’il convenu d’appeler l’estimateur de Moindre Carrés Quasi Générali-
sés (MCQG). L’estimateur de MCQG est similaire a celui de MCG a ceci puisque l’on remplace la
matrice A inconnue par un estimation convergente
63
2
σ̂w
Avec θ̂ = (σ̂2 +T σ̂v2 )
La mise en ouvre de cette méthode nécessite donc de disposer l’estimation de
w
variances inconnues σu2 etσw2 . Il existe différente façon d’estimer les variances σu2 etσw2 . Selon Swamy et
ε̂0w εw ˆ0 w εw
Arora (1972), l’estimateur naturel de la varianceσw2 est donnée par σ̂w2 = rang(M w)
= N (Tε−1)−K w
Effects :
var std.dev share
idiosyncratic 0.003613 0.060105 0.202
individual 0.014296 0.119567 0.798
theta : 0.8713
Residuals :
Min. 1st Qu. Median 3rd Qu. Max.
−0.1380 − 0.0391 − 0.0052 0.0372 0.1780
Coefficients :
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 9.622504 0.208458 46.1604 < 2.2e − 16 ∗ ∗∗
log(output) 0.905881 0.025391 35.6769 < 2.2e − 16 ∗ ∗∗
log(pf) 0.423126 0.013981 30.2640 < 2.2e − 16 ∗ ∗∗
lf −1.064454 0.200456 − 5.3102 8.459e − 07 ∗ ∗∗
— Signif. codes : 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
64
c) Tests de l’existence de spécification individuelle
*Test de Brench-Pagan : Le test de Brench-Pagan a pour but de tester l’absence des effets
spécifiques individuels. Il s’agit d’un test de multiplicateur de Lagrange qui consiste à tester la nullité
de leur variance :
(
H0 : σu2 = 0
H1 : σu2 6= 0
Sous l’hypothèse nulle, il n’y a pas d’effet aléatoire individuel
*Test d’Hausman : Dans le modèle avec les effets individuels. Les effets individuels peuvent
être soit fixes ou aléatoires. Le test Hausman étudier la présence d’une corrélation entre les effets
spécifiques et les variables explicatives. Sous l’hypothèse nulle les régresseurs sont strictement exo-
gènes, les estimateurs de MCO et de MCG sont alors convergent mais l’estimateur de MCQG est
non convergent. Dans ce cas, on doit appliquer l’estimateur de MCO pour estimer le modèle à effets
fixes individuels.
Soient β̂M CG et β̂M CQG les estimateurs deβ. Test d’Haussmann consiste à calculer le statistique de
Wald
W = (β̂M CG − β̂M CQG )0 vâr(β̂M CG − β̂M CQG )−1 (β̂M CG − β̂M CQG )
Avec
vâr(β̂M CG − β̂M CQG ) = vâr(β̂M CG ) − vâr(β̂M CQG )
65
Cette statistique est distribuée selon une loi du χ2 de k degré de liberté. Une valeur de W supé-
rieur à la valeur de χ2 (k)conduit à rejeter l’hypothèse nulle et à considère l’estimateur de MCQG est
convergent (le modèle à effets aléatoires).
phtest(air.fe,air.re)
Hausman Test
data : log(cost) log(output) + log(pf) + lf
chisq = 2.7629, df = 3, p-value = 0.4296
alternative hypothesis : one model is inconsistent
Test d’Hausman consiste à retenir ou pas l’hypothèse nulle Ho : Il y a corrélation entre effet spécifique
et les variables explicatives. Règle de décision : Si W < χ2 (k) on accepte l’hypothèse nulle. D’après
le résultat on a statistique Wald W = 2.7629; > χ2 (3) = 0.35 A un niveau de significativité 5% ; on
rejette H0 selon la quelle les effets individuels sont corrélés avec les variables explicatives. Par suite
le modèle étudié est à effets aléatoires individuels.
Exercice Nous nous intéressons à l’évolution des dépenses d’éducation sur 25 ans (variable à
expliquer notée y) pour neuf pays. Les facteurs explicatifs proposés sont :
x1it = évolution des dépenses militaires en milliers de dollars pour le pays i à l’année t,
x2it = évolution du PIB (Produit intérieur brut) en milliers de dollars pour le pays i à l’année t.
Le modèle à étudier est donc :
yit = a0i + a1i x1it + a2i x2it + it pouri = 1, ..., N ett = 1, ..., T
66
Chapitre 7
La régression linéaire (simple ou multiples) est utilisée lorsque la variable dépendant (Y) est
continue (échelle d’intervalle ou du rapport) comme par exemple pour estimer le prix d’une résidence
en fonction de la superficie, le nombre de salles de bain dans une résidence. Cependant, la régression
logistique est utilisée lorsque la variable dépendante (Y) est catégorielle (échelle nominales), comme
par exemple pour prédire la réussite à un cours par des variables dépendantes qui peuvent être
continues ou catégorielles.
Le modèle de la régression logistique est utilisé dans des nombreux domaines. En médecine, elle
permettra par exemple de trouver les facteurs qui caractérisent un groupe de sujets malades par
rapport à des sujets sains. Dans le domaine des assurances, elle permettra de cibler une fraction de la
clientèle qui sera sensible à une police d’assurance sur tel ou tel risque particulier. Dans le domaine
bancaire, elle permettra de détecter les groupes à risque lors de la souscription d’un crédit.
L’objectif de la régression logistique est le même que celui de la régression linéaire ; il s’agit
d’analyser une variable dépendante catégorielle (échelle nominales) Y avec une variable dépendante
X. En d’autre termes, la régression permet de déterminer les effets des variables X sur la variable
dépendante Y. Cette analyse permet de mieux comprendre et d’expliquer les phénomènes étudiés,
voir même de les prédire, si le modèle est suffisamment solide.
67
7.1.1 Formulation mathématique
Pour chacune d’expérience, la variable endogène observée yi prend une valeur 0 ou 1 selon la
réalisation d’un événement (l’individu a ou non bien supporté l’expérience). Il est clair qu’on ne peut
pas utiliser la loi normale dans la modélisation de ce type des variables qualitatives. En fait, comme
yi ne peut prendre que deux valeurs (0 ou 1), il en est de même de la perturbation εi prend la valeur
1 − βxi avec une probabilité pi ou la valeur −βxi avec une probabilité pi .
Par suite, la perturbation εi admet obligatoirement une loi discrète, ce qui interdit l’hypothèse
de la normalité. Par conséquent, le modèle linéaire classique n’est pas adéquat pour formaliser la
dépendance de la variable endogène de nature discrète vis-à-vis des valeurs prises par des facteurs
explicatifs. La solution proposée est de modéliser la moyenne conditionnelle de la variable à expliquer
étant donnée les variables explicatives µi = E(yi /xi ) au lieu de la variable à expliquer yi elle-même.
Comme les données de la variable à expliquer yi sont dichotomiques, alors Y suit la loi de Bernoulli
de paramètre µi = P (yi = 1). Pour modéliser cette probabilité on suppose que la décision repose sur
la valeur prise par une variable inobservable yi∗ appelée variable latente, selon le schéma suivant :
(
On observe yi = 1 lorsque yi∗ ≥ 0
On observe yi = 0 lorsque yi∗ < 0
En réalité, on ne dispose pas des informations sur la variable latente yi∗ qui permet à l’individu
de prendre la décision (choix 1 ou 0). Pour rendre le modèle estimable, on suppose que cette variable
latente dépend linéairement d’un certain nombre de variables explicatives.
yi∗ = xi β + εi
Les perturbations εi sont supposées indépendantes, de moyenne nulle et suivent une même loi
ayant une fonction de répartition F. Cette hypothèse d’indépendance se traduit par la condition que
les observations doivent être faites sur des individus différents.
On déduit facilement le lien entre le prédicteur linéaire βxi et la moyenne µi = E(yi /xi ).
On peut alors écrire : µi = E(yi /xi ) = P (yi = 1/xi ) = P (yi ≥ 0/xi ) = P (xi β + εi ≥ 0/xi ) = P (εi ≤
xi β/xi ) = F (xi β)
La fonction de répartition F dépend alors de l’hypothèse faite sur la distribution des perturbations
εi . On retient habituellement pour cette distribution soit une loi normale centrée réduite (le modèle
est appelé modèle Probit), soit une loi logistique ayant une distribution centrée et de variance (π/3)(le
modèle est appelé modèle Logit). La fonction de répartition associée à la loi logistique s’écrit sous la
forme suivante :
1
µi = Λ(xi β) =
1 + exp −(xi β)
On définie la fonction logit g : [0, 1] → R est une fonction réciproque de la fonction de répartition
de la loi logistique Λ(xi β) tel que :
µi
g(µi ) = log = xi β
1 − µi
68
Le modèle logistique n’admet pas une variable expliquée (comme dans le cas de la spécification
linéaire) mais le codage quantitatif associés à la réalisation d’un événement.
N
Y N
Y
1−yi
L(y, β ) = µyi i (1 − µi ) = [F (xi β)]yi [1 − F (xi β)]1−yi
i=1 i=1
Il ne reste plus alors qu’à spécifier la fonction de distribution F (.) pour obtenir la forme fonc-
tionnelle de la vraisemblance.
Ainsi, dans le cas du modèle logit, ∀xi β ∈ R , on a :
e xi β
F (xi β) = .
1 + exβi
A partir de cette définition, on déduit alors la log -vraisemblance comme suit :
N
X
log L(y, β) = yi log[F (xi β)] + (1 − yi ) log(1 − F [xi β)]
i=1
En distinguant les observations yi = 1et celles pour les quelles on a yi = 0, le log -vraisemblance
peut s’écrit sous la forme :
X X
Log L(y, β) = log F (xi β] + log(1 − F (xi β))
i:yi=1 i:yi =0
N
∂ log L(y, β) X f (xi β) 0 f (xi β) 0
G(β) = = yi xi + (yi − 1) xi
∂β i=1
F (xi β) 1 − F (xi β)
69
0
Où f (.) est la fonction de densité associé à F (.) et où xi désigne le transposée du vecteur xi de
dimension (1 ; k). En simplifiant, l’expression du gradient, on obtient alors :
N
X yi − F (xi β) f (xi β) 0
G(β) = x
i=1
F (xi β) (1 − F (xi β)) i
On peut en outre exprimer le gradient en distinguant les observations yi=1 et celles pour lesquelles
on a yi=0
X f (xi β) 0 X f (xi β) 0
G(β) = xi − xi
i:y
F (xi β) i:y =0
1 − F (xi β)
i=1 i
70
7.1.4 Propriétés asymptotiques des estimateurs de maximum de vrai-
semblance
Lorsque l’on cherche à établir les propriétés asymptotiques des estimateurs du maximum de
vraisemblance dans le cadre de modèles logistiques, et plus généralement dans le cadre de modèle
à variables qualitatives, toute la difficulté réside dans le fait que l’on ne dispose pas d’expression
analytique pour ces estimateurs.
En effet, nous avons vu que les équations de vraisemblance associée au logit sont non linéaires
dans les paramètres. Dés lors, il n’est pas possible alors d’exprimer les estimateurs, solutions de ces
équations, comme des fonctions simples des observations. Nous avons vu q’il était alors nécessaire de
recouvrir à des algorithmes d’optimisation numérique.
Mais devant l’impossibilité d’écrire les estimateurs du maximum de vraisemblance comme des
fonctions simples d’observations, il est alors difficile d’étudier la convergence de ces estimateurs
comme nous avions pu le faire dans le cas des modèles linaires standard. Il convient ainsi d’adopter
une démarche particulier ou l’on va chercher à étudier la convergence du critère de maximum de
vraisemblance, afin de montrer la convergence des estimateurs du MV, solution du programme de
maximisation de ce critère.
Un certain nombre de rappels sur les différentes notions de convergence sont proposés Toute-
fois, la lecture de ces rappels doit nécessairement s’accompagner d’une étude plus systématique des
fondements probabilistes de ces notions.
a)Convergence du critère de MV
b) Qualité de l’estimateur
71
la significativité du modèle c’est à dire l’hypothèse H0 : β2 = β3 = ... = βk = 0 dans un modèle
comportant k paramètres on peut appliquer le test du ratio vraisemblance.
∗
LR = −2(log L(β̂M V ) − logL(β̂M V ))
∗
P
Où log L(β̂M V) = yi log F (βb1 ) + (1 − yi )log (1 − F (βb1 ))
et logL(β̂M V ) est la valeur de la fonction log L lorsque les paramètres βk = 0..
Sous l’hypothèse H0 , pi = F (β1 ) = P la probabilité pi est identique pour tous les individus et elle
est donnée par la proportion P d’individus pour lesquels on observe yi = 1
Si la statistique LR est inférieure au χ2α (k − 1) on accepte l’hypothèse H0 .
Par ailleurs, pour mesure la quantité de l’agissement, il existe plusieurs possibilités une manière
simple consiste à calculer le nombre de mauvaise prédictions. L’estimateur du maximum de vraisem-
blance donne p̂i . On suppose que si p̂i est inférieure à 0,5 alors la valeur estimée de yi notée ŷi est
égale à 0 si non elle est égale à 1.
On peut également utiliser le coefficient de détermination (R2 ) suggéré par MC Fadden RM 2
F
2 log L(β̂M V )
RM F = 1− ∗
logL(β̂M V)
L’estimation du modèle logit permet d’évaluer la probabilité Pi et l’influence des différentes
variables explicatives sur cette probabilité.
Dans le modèle linéaire les paramètres estimés mesurent les effets marginaux c’est à dire l’impact
de la variation des variables explicatives sur la variable endogène. Dans les modèles logit les effets
marginaux ne sont pas égaux aux valeurs estimées des paramètres.
0 L
[g(β̂) − r] [GV̂ (β̂)G0 ]−1 [g(β̂) − r] −→ χ2 (c) si N → ∞
AvecG = ∂g(.)
∂β
et V̂ (β̂ ) l’estimateur de la matrice de variance covariance des coefficients. Dans le
cas qui nous intéresse, on a g(β) = βj et r = a et vecteur G, de dimension (K, 1) comporte k-1
72
zéros et 1 à la j me position. La statistique du test de Wald associé au test unidirectionnelH0 : βj = a
contre H1 : βj 6= a admet la loi suivante sous H0 :
2
h i0 h i β̂j − a L
β̂ − a (v̂jj )−1 β̂j − a) = −→ χ2 (1)
V̂ii
Où v̂jj désigne l’estimateur de la variance de l’ j coefficient βj . Ainsi si l’on note χ295% (1) la
me
(β̂−a) 2
quantité à 95% de la loi Khi deux. Le test de Wald au seuil 5% consiste à accepter H0 si ν̂jj est
inférieur χ295% (1) et à refuser H0 si cette quantité à supérieur àχ295% (1).
La plupart des logiciels ne proposent pas cette statistique de Wald mais une statistique Zj définie
comme la racine carrée de la précédente. Compte tenue du lien entre la loi normal centrée réduite et
la loi du Khi deux à un degrés de liberté, on a immédiatement sousH0 .
β̂j − a L
Zj = p −→ N (0, 1)
Vjj
Test du rapport des maxima de vraisemblance Dans le cas des modèles dichotomiques, on
peut appliquer sans difficulté particulière la logique du test du rapport du maximum de vraisemblance.
Ainsi, on estime le modèle non contraint et d’autre part le modèle contraint : soient βj1 et βjc les deux
estimateurs. Ainsi, on obtenue la statistique LRT en calculant l’écart de la log vraisemblance. La
statistique LRTj du test du rapport des maxima de vraisemblance associé au test unidirectionnel
H0 : βj = a contre H0 : βj 6= a admet la loi suivante sous H0 :
L
LRTj = −2 [ log L(y, β̂j) − log L(y, βbjc )] −→ χ2 (1)
Où β̂j et β̂jc désignent respectivement les estimateurs non contraint et contraint de βj .
Le test du rapport des maxima de vraisemblance au seuil de 5% consiste à accepter H0 si LRTj <
χ95% (1)et à refuser H0 si LRTj ≥ χ295% (1). Cette procédure est asymptotiquement équivalente à celle
2
7.2.1 Les effets marginaux lorsque les variables explicatives sont continues
Pour mesurer l’impact d’une variation d’une variable explicative sur la probabilité on doit calculer
δPi 0
= βf (xi β)
δxi
73
Si le modèle estimé est un modèle probit
δ P̂i 0
= β̂probit φ(xi β̂probit )
δxi
Si c’est un modèle logit
δ P̂i 0 0
= β̂logit Λ(xi β̂logit )(1 − Λ(xi β̂logit ))
δxi
1
Avec Λ = 1+exp(−x) Ces effets marginaux peuvent être calculés pour chaque individu. Cependant pour
simplifier la présentation des résultats, on évalue souvent ses effets marginaux au point moyen x̄ ou
pour un niveau de probabilité initial donné.
7.2.2 Les effets marginaux lorsque les variables explicatives sont qualita-
tives
Si dans un modèle logit il ya une variable explicative qualitative, on ne peut pas appliquer la
procédure précédente pour mesurer l’impact sur la probabilité Pi . Par exemple, si on fait intervenir
la variable explicative x2 indiquant le sexe de l’individu x2 = 1 si l’individu i est de sexe masculin
x2 = 0 si l’individu i est de sexe féminin .
Le modèle s’ecrit alors
Pi = F (β0 + β1 x1i + β2 x2i )
Et son estimation avec la méthode du maximum de vraisemblance donne .
Au point moyen x¯1 = x1i pour évaluer l’impact de la variable x2 sur la probabilité d’appartenir à
la classe 1, on calcul la probabilité q’un individue de sexe masculin on la note PˆM et pour un individu
de sexe féminin on la note PˆF . Dans le cas d’un modèle logit
7.3 Application
Dans cet exemple, nous allons ajusté un modèle de régression logistique à un ensemble de don-
nées relatives à la santé. Nous avons des données représentant 8 personnes atteintes de cancer et 8
personnes en bonne santé et le nombre de cigarettes fumées par chacun d’eux et son sexe.
La variable Diagnostic est une variable binaire qui prend la valeur 1 pour les individus qui ont
atteintes de cancer, et la valeur 0 en bonne santé.
Le modèle de régression logistique permet de déterminer l’effet de nombre des cigarettes fumés et la
74
nature selon le sexe sur la probabilité d’être atteint du cancer. Nous allons utiliser la commande glm
() pour exécuter une régression logistique, en régressant le Diagnostic sur les nombre de cigarette
et sur le gendre.
Diagnostic=c(0,0,0,0,1,0,0,0,1,1,1,1,1,0,1,1)
sexe=c(2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1)
noCig =c(10,13,12,12,7,8,12,17,32,37,28,34,32,37,30,28)
myData=data.frame("noCig"=noCig,"sexe"=sexe,"Diagnostic"=Diagnostic)
mod1=glm(cancer ∼ noCig+gender,data=myData,family=binomial)
summary (mod1)
exp(cbind("OR"=coef(mod1),confint(mod1)))
Dans la formule entre les parenthèses, à gauche de ∼ se trouve la variable dépendante Diagnostic
qui doit être codé 0 et 1 pour que glm le lise en tant que binaire. A droite de ∼, nous listons les
deux variables prédictives. Enfin, après la virgule, nous précisons que la distribution est binomiale.
La fonction de liaison par défaut dans glm pour une variable de résultat binomiale est le logit. Nous
pouvons accéder à la sortie du modèle en utilisant summary (mod1).
Le modèle estimé génère la relation suivante entre le logit (log odds) et les deux prédicteurs :
p
logit( ) = −7.8 + 0, 18 ∗ noCig + 2.55 ∗ sexe
1−p
Le coefficient de noCig est égal à 0.18, de sorte qu’une unité de cigarette supplémentaire entraîne
une modification d’environ 0,18 unité la cote d’être atteint de cancer (c’est-à-dire un changement
de 0.18 unité dans le logit). Les signes des deux prédicteurs montrent que noCig influe de manière
positive sur l’atteint de cancer. En outre, nous pouvons convertir ces logits en rapports de cotes.
Nous constatons que OR( noCig) = 1,19798 et donc l’augmentation du nombre de cigarettes
conduit à une augmentation de la probabilité d’être atteint de cancer de 19,79% et avec un intervalle
de confiance à 95% pour le rapport de cotes des cigarettes sont [1.0512,1.5269].
La sortie produite par glm () inclut plusieurs quantités supplémentaires qui nécessitent une dis-
cussion. Nous examinons la valeur z pour chaque estimation. La valeur z est la statistique de Wald qui
teste l’hypothèse selon laquelle l’estimation est égale à zéro. L’hypothèse nulle est que l’estimation a
une distribution normale avec un zéro moyen et un écart type de 1. La valeur p indiquée, P (> |z|),
donne la zone d’extrémité dans un test bilatéral. Nous constatons que le nombre de cigarettes avait
un effet significatif sur l’incidence du cancer
Pour notre exemple, nous avons une déviance nulle d’environ 22.181 sur 15 degrés de liberté.
Cette valeur indique un mauvais ajustement (différence significative entre les valeurs ajustées et les
75
valeurs observées). L’inclusion des variables indépendantes (Nocig et sexe) a réduit la déviance de
près de 9.4 points sur 2 degrés de liberté. La déviance résiduelle est de 12.772 sur 13 degrés de liberté
(c’est-à-dire une perte de deux degrés de liberté).
Nous allons maintenant créer un graphique pour chaque prédicteur. Cela peut être très utile
pour nous aider à comprendre l’effet de chaque prédicteur sur la probabilité d’une réponse de 1 sur
notre variable dépendante. Nous souhaitons tracer chaque prédicteur séparément, nous ajustons donc
d’abord un modèle distinct pour chaque prédicteur. Ce n’est pas la seule façon de le faire, mais un
moyen que j’estime particulièrement utile pour décider quelles variables doivent être entrées comme
variables prédites.
Maintenant, nous créons nos parcelles. Nous établissons d’abord une séquence de valeurs de
longueur que nous utiliserons pour tracer le modèle ajusté. Une séquence de 0 à 15 convient à
merveille pour calculer les calculs. Le modèle a produit une courbe indiquant la probabilité que
diagnostic = 1 pour le nombre de cigarette. Clairement, plus le nombre de cigarette est élevé, plus
la personne risque d’atteindre le cancer.
Notez que la seule variable qui mène au cancer est le nombre de cigarettes p = 0,0313 <0,05.
L’introduction de la variable du nombre de cigarettes a montré une amélioration du modèle car la
déviance résiduelle était inférieure à Null Deviance et AIC, bien que petite, ne pourront pas utiliser
pour l’évaluation du modèle, car la norme AIC n’apparaît utile que lorsque plusieurs modèles sont
comparés.
76
Chapitre 8
La régression de Poisson permet de modéliser des comptages distribués selon une loi de Poisson en
fonction de variables explicatives quantitatives ou qualitatives. La régression de Poisson est similaire
à la régression multiple régulière sauf que la variable dépendante (Y) est un compte observé qui suit
la distribution de Poisson. Ainsi, les valeurs possibles de Y sont les entiers non négatifs : 0, 1, 2, 3, etc.
On suppose que les grandes quantités sont rares. Par conséquent, la régression de Poisson est similaire
à la régression logistique, qui comporte également une variable de réponse discrète. Cependant, la
réponse ne se limite pas à des valeurs spécifiques comme dans la régression logistique.
Un exemple d’application appropriée de la régression de Poisson est une étude de la relation
entre le nombre de colonies de bactéries et les diverses conditions environnementales et dilutions.
Un autre exemple est le nombre de pannes pour une machine donnée dans différentes conditions de
fonctionnement. Un autre exemple encore est celui des statistiques vitales concernant la mortalité
infantile ou l’incidence du cancer parmi des groupes de différentes démographies.
e−µ µy
P (Y = y/µ) = (y = 0, 1, 2, ..)
y!
Notez que la distribution de Poisson est spécifiée avec un seul paramètre µ. C’est le taux d’incidence
moyen d’un événement rare par unité d’exposition. L’exposition peut être le temps, l’espace, la
distance, la superficie, le volume ou la taille de la population. Comme l’exposition est souvent une
période de temps, nous utilisons le symbole t pour représenter l’exposition. Lorsqu’aucune valeur
d’exposition n’est donnée, on suppose qu’il en est une.
Le paramètre µ peut être interprété comme le risque d’une nouvelle occurrence de l’événement
pendant une période d’exposition spécifiée, t. La distribution de Poisson a la propriété que sa moyenne
et sa variance sont égales.
77
8.2 Le modèle de régression de Poisson
Dans la régression de Poisson, on suppose que le taux d’incidence de Poisson µ est déterminé par
un ensemble de k variables de régression (les X). L’expression relative à ces quantités est
Les équations de vraisemblance peuvent être formées en prenant les dérivées par rapport à chaque
coefficient de régression et en fixant le résultat à zéro. Cela conduit à un ensemble d’équations non
linéaires qui n’admet aucune solution de forme fermée.
Les équations de vraisemblance sont :
n
∂lnL X
= (yi − µi )xi = 0
∂β i=1
Le Hessien est défini négatif pour tout et tout x et β. La méthode de Newton est alors pour ce modèle
un algorithme simple qui converge rapidement.
78
Ainsi, un algorithme itératif doit être utilisé pour trouver l’ensemble des coefficients de régression
maximisant la log-vraisemblance. En utilisant la méthode des moindres carrés pondérés de manière
itérative, une solution peut être trouvée en cinq ou six itérations. Cependant, l’algorithme nécessite
un passage complet des données à chaque itération, il est donc relativement lent pour les problèmes
avec un grand nombre de lignes. Avec les ordinateurs d’aujourd’hui, cela devient de moins en moins
problématique.
On en déduit la matrice de variance-covariance asymptotique estimée de l’estimateur du maximum
0 0
de vraisemblance :−[ ni=1 −µ̂i xi xi ]−1 La prédiction pour l’observation i est µ̂i = exp(xi β̂)
P
Rappelez-vous que dans le modèle de Poisson, la moyenne et la variance sont égales. En pratique, les
données rejettent presque toujours cette restriction. Habituellement, la variance est supérieure à la
moyenne - une situation appelée surdispersion. L’augmentation de la variance est représentée dans
le modèle par un multiple constant de la matrice de variance-covariance.
Ces deux statistiques sont approximativement distribuées en khi carré avec n - k degrés de liberté.
Lorsqu’un test est rejeté, il y a un manque d’ajustement important. Lorsqu’un test n’est pas rejeté,
il n’y a aucune preuve de manque d’ajustement. La statistique Pearson est uniquement distribuée
en khi-carré lorsque vous analysez des données groupées. Par conséquent, si vous n’utilisez pas de
variable de fréquence, vous ne devez pas utiliser la statistique de Pearson comme test d’ajustement.
La statistique de Pearson est souvent utilisée comme test de surdispersion.
8.3.3 Déviance
La déviance est le double de la différence entre le log-vraisemblance maximum réalisable et le
log-vraisemblance du modèle ajusté. Dans la régression multiple sous normalité, la déviance est la
79
somme résiduelle des carrés. Dans le cas de Régression de Poisson, la déviance est une généralisation
de la somme des carrés. La formule pour la déviance est
LLf it − LL0
R2 =
LLmax − LL0
La statistique R-carré ne s’étend pas aux modèles de régression de Poisson. Différents tests pseudo-
R-carrés ont été proposés. Ces pseudo-mesures ont pour propriété, lorsqu’elles sont appliquées au
modèle linéaire, de correspondre à l’interprétation du modèle linéaire R au carré. Dans la régression
de Poisson, la mesure pseudo-R au carré la plus populaire est fonction des log-vraisemblances de
trois modèles.
80
8.5 Application
Dans cet exemple, nous ajustons un modèle de régression de poisson pour modéliser les données
des comptages. L’ensemble de données comprend Y le nombre d’élèves du secondaire ayant atteint
une maladie infectieuse au cours d’une période de plusieurs jours à compter du début de l’épidémie.
Comme Y est une variable de comptage à valeurs entières, on peut supposer que Y ∼ P (λ(x)),
où λ(x) désigne le nombre moyen inconnu d’infections. On veut estimer λ(x) à partir des données.
On adopte le modèle de régression de Poisson (avec interactions si on travaille avec des variables
qualitatives) :
ln(λ(x)) = β0 + β1 Days
Objectifs : Estimer les paramètres inconnus à partir des données et étudier la qualité du modèle.
La modélisation de la régression de Poisson et les estimations des paramètres par la méthode du
maximu de vraisemblance s’obtiennent par les commandes de R suivante :
model1 =glm(Students ∼ Days, poisson)
summary(model1)
On a les estimations ponctuelles des coefficients β0 , β1 , ainsi que les degrés de significativité des va-
riables considérées sur Y .
Le coefficient (β1 = −0.0174) est négatif ce qui indique qu’au fur et à mesure que le nombre de jours
de diagnostic augmente, le nombre moyen d’élèves atteints de la maladie diminue. Ce coefficient est
hautement significatif (p <2e-16).
La commande predict() permet d’estimer le nombre étudiant ayant atteint la maladie pour le
10 jours. La valeur prédite de Y quand Days = 10 :
predict.glm(model1, data.frame(Days=10), type = "response")
81
Ainsi, le nombre moyen d’infections pour un individu vérifiant les modalités ( Daus=10) est de
6.14. La probabilité que nombre de étudiant ayant atteint la maladie est égale à 5 est donnée par :
probs = dpois(5, 6.14))
probs
Cela renvoie : p-valeur= 1 > 0.05, donc le modèle est bien adapté au problème.
Ceci est une hypothèse de modèle très importante. Dans la section suivante, nous allons donc adapter
le modèle à l’aide d’erreurs de type quasi poisson.
82
√
inchangés, mais leurs écarts-types seront multiplés par φ. Ainsi, certains paramètres qui étaient
marginalement significatifs peuvent ne plus le rester.
La dispersion excessive pose problème si la variance conditionnelle (variance résiduelle) est supé-
rieure à la moyenne conditionnelle. Un moyen de vérifier et de traiter la dispersion excessive consiste
à utiliser un modèle de quasi-poisson, qui intègre un paramètre de dispersion supplémentaire pour
prendre en compte cette variance supplémentaire.
Nous adaptons un modèle quasi-Poisson aux mêmes données. Dans ce cas, il est d’usage d’utiliser
family=« quasipoisson » plutôt que « poisson ». La différence est que avec quasi poisson le paramètre
de dispersion n’est plus fixé à 1 (dans l’exemple, il va être de 6.3 environ ).
Dans R, la famille ’quasipoisson’ peut être utilisée pour traiter ces problèmes de surdispersion (de
la même manière la famille quasibinomial’ peut être utilisée). L’estimation de φ sera donné dans le
résumé du modèle GLM quasi Poisson. En plus, les coefficients estimés sont inchangés, mais leurs
écart-type sont corrigés (et avec les statiques de test et p-value associées) pour tenir compte de ce
phénomène d’overdispersion, comme le montrent les résultats ci-dessous.
Le résultat de notre tentative de rendre compte de la dispersion excessive est que la déviance
résiduelle n’a pas changé. Le paramètre de dispersion, qui devait obligatoirement être 1 dans notre
dernier modèle, peut être estimé ici. En fait, il est estimé à 0,79. Ce paramètre nous indique combien
de fois la variance est plus grande que la moyenne. Comme notre dispersion était inférieure à un,
il s’avère que la variance conditionnelle est en réalité inférieure à la moyenne conditionnelle. Nous
avons une sous-dispersion, pas terminée.
Quoi qu’il en soit, nous traçons maintenant la régression. Nous établissons un axe temporel allant
de 0 à 150 (le nombre de jours). Cependant, nous incluons de petites augmentations de 0,1 afin de
donner une apparence lisse à notre graphique. Nous allons évaluer le modèle sur ces valeurs, puis
utiliser ces valeurs pour tracer le modèle.
83
timeaxis =seq (0,150,0.1)
Y =predict(model2, list(Days = timeaxis)))
plot(Days, Students, xlab = "DAYS", ylab = "STUDENTS", pch = 16)
lines(timeaxis, exp(Y), lwd = 2, col = "blue")
Bien sûr, au lieu de prendre l’exponentielle des valeurs ajustées, nous aurions aussi pu utiliser
la fonction Predict () avec l’argument type = "réponse". Calculons l’impact sur le nombre de cas
résultant d’une augmentation d’un jour sur l’axe des temps.
Enfin, nous traçons le modèle ajusté et nous prenons l’exponentielle des valeurs ajustées car les
valeurs ajustées sont renvoyées sur une échelle logarithmique. Prendre les rétro-transformations ex-
ponentielles de l’échelle du log des données d’origine.
Le graphique montre une diminution non linéaire du nombre de cas avec le nombre de jours.
Nous pouvons calculer l’évolution du nombre d’élèves présentant la maladie pour chaque jour
supplémentaire, comme suit :1 - 0,9826884=0,0173116
La réduction (rapport de taux) est d’environ 0,02 cas par jour supplémentaire.
En examinant le résumé du modèle nous pouvons voir que φ est estimé à 0.793. Nous avons donc eu
raison de corriger le modèle afin de prendre en compte la surdispersion. Par contre, si nous regardons
la significativité du coefficient de regression associé à l’élévation, nous remarquons qu’il n’est plus
significatif. Cependant, 0.793 indique que la surdispersion est forte et en général un GLM quasi-
Poisson est favorisé lorsque φ est compris entre 1 et 15. Dans un soucis pédagogique, nous allons
considérer que nous ne sommes pas satisfait avec ce modèle et que nous voulons maintenant ajuster
une distribution binomiale négative à nos données. Deux points sont importants à garder en tête
lorsque vous utilisez un GLM quasi-Poisson afin de corriger la surdispersion :
-Les GLMs quasi-Poisson n’ont pas d’AIC. En effet, la vraisemblance d’un modèle GLM quasi-Poisson
84
ne peut pas être spécifiée et s’appuie sur une procédure de pseudo-maximum de vraisemblance. Par
conséquence les GLMs quasi-Poisson n’ont pas d’AIC, et ce critère ne peut pas être utilisé afin de
comparer différents modèles. Toutefois des alternatives ont été developpées pour gérer cette situation
(e.g. quasi-AIC). -La surdispersion influence la comparaison de modèles. En effet, la surdispersion
influence la comparaison de deux modèles emboités et doit donc être prise en considération. Par
exemple, considérons que nous voulons comparer le modèle GLM1, qui contient p1 paramètres avec
le modèle GLM2, qui contient p2 paramètres. GLM1 est emboité dans GLM2 et p2 > p1. La compa-
raison des deux modèles est basées sur le test du rapport des vraisemblances des deux modèles, D1
et D2 respectivement. Si la surdispersion est connue, les déviances doivent être corrigées de manière
approprié selon D∗ = Dφ , et le test final sera basé sur le critère D1* - D*2 qui est supposé être dis-
tributé selon une distribution du avec p2-p1 degrés de liberté lorsque le modèle GLM1 est correct.
Mais dans certain cas φ n’est pas connu. Par exemple, lorsque vous spécifiez un GLM avec un distri-
bution normale. Dans ce cas, φ peut être estimé a posteriori en utilisant la déviance résiduelle du plus
gros modèle de telle sorte que le critière de comparaison devienne [(D1-D2)/(p2-p1)]/[D2/(n-p2)].
Ce critère est supposé suivre une distribution F avec p2-p1 et n-p2 degrés de liberté.
De cette manière nous pouvons voir comment cette distribution va gérer la surdispersion dans
les modèles GLM. Le deuxième terme de la variance de la distribution BN va déterminer le degré
de surdispersion. En effet, la surdispersion est indirectement déterminée par k, que représente le
2
paramètre de dispersion. Si k est grand (par rapport à µ2 ), la deuxième partie de la variance, µk
va s’approcher de 0, et la variance de Y sera µ. Dans ce cas la distribution BN converge vers la
distribution de Poisson et vous pourriez tout aussi bien utiliser cette dernière. Par contre, plus k sera
petit et plus la surdispersion sera grande. Comme avec toutes les autres distributions, un GLM avec
une distribution BN se spécifie en trois étapes. Tout d’abord le modèle fait l’hypothèse que les Yi
suivent une distribution BN de moyenne µi et de paramètre k.
85
Les deux dernières étapes définissent le prédicteur linéaire ainsi que la fonction de lien entre la
moyenne des Yi et le prédicteur linéaire. La fonction de lien utilisée par les GLMs avec une distribution
BN est le logarithme ce qui permet de s’assurer que les valeurs prédites soient toujours positives.
La distribution NB n’est pas dans la fonction glm(), il faut donc installer et charger la librairie
MASS. Vous pouvez ajuster un GLM avec une distribution BN à l’aide de la fonction glm.nb() :
library("MASS")
glm.negbin = glm.nb(Students ∼ Days)
lsummary(glm.negbin)
Le résumé du modèle et similaire à celui des autres GLMs (e.g. GLMs Poisson). Cependant
vous avez maintenant un nouveau paramètre, theta, qui est le paramètre k de la variance de votre
distribution. L’écart-type de ce paramètre est aussi fourni, mais attention à son interprétation car
l’intervalle n’est pas symétrique.
Tous les GLMs que nous venons de voir pour modéliser des données de d’abondance (Poisson, quasi-
Poisson et BN) utilisent la même relation log-linéaire entre moyenne et prédicteur linéaire (log(µ) =
X.β). Toutefois ils vont autoriser différentes relations entre la moyenne et la variance et vont aussi se
reposer sur des méthodes d’estimation de la vraisemblance différentes. Les GLMs Quasi-Poisson ou
BN sont privilégiés afin de traiter la surdispersion. Malheureusement dans certains cas les données
peuvent contenir trop de zeros et d’autres modèles seront plus éfficaces pour traiter ces situations.
C’est par exemple le cas des “zero-augmented models” (e.g. zero-inflated Poisson [ZIP]) qui vont
traiter les zéros indépendamment des autres valeurs.
86
Chapitre 9
87
théorème de Bayes, ainsi que les principales de distributions des probabilités à priori.
Soit θ une matrice des paramètres a estimé à partir d’un échantillon d’observations y, alors, d’après
la loi de Bayes, nous avons :
p(y|θ)p(θ)
p(θ|y) = p(y)
,
p(θ|y) v (θ)p(y|θ).
On utilise l’équivalence, pour éviter la normalisation à partir de la quantité p(y). Il faut noter que
la quantité de vraisemblance contient l’information de l’échantillon à partir des observations générées
dans le temps.
Toutefois, la distribution à priori présente une loi supposée par le modélisateur afin de déterminer
la distribution à posteriori après une combinaison avec l’information disponible.
Plusieurs articles de recherche sont impliqués sur la question de ces distributions à priori, ainsi
qu’il existe plusieurs travaux récents sur le thème. Toutefois, notre objectif au niveau de cette partie,
est de pouvoir introduire les notions de base sur les distributions à priori.
On a choisi à partir de nos recherches théoriques de présenter les distributions à priori les plus
connues, afin de pouvoir nous introduire à cette notion, qui sera utile pour nos travaux antérieurs.
88
La distribution a priori uniforme
C’est la première distribution à priori supposée par Bayes. En effet, il considère que l’absence de
l’information sur le comportement du paramètre, nous pousse à supposer une loi à priori uniforme
pour ce dernier.
Cette approche était fortement critiquée par les auteurs tels que, Fisher, Keynes. . . qui considèrent
que cet argument d’absence de l’information est critiquable.
L’idée dans cette approche est de pouvoir mettre une relation entre la distribution à priori du
modèle et la quantité d’informations de Fisher de ce dernier. En effet, notons par f la vraisemblance
de x|θ avec θ ∈ Rd . L’information de Fisher est donnée par :
L’idée de Jeffrey est de construire une distribution à priori uniforme en se basant sur le critère
de l’information de Fisher. Par conséquent, la distribution à priori donnée par cette approche est la
suivante :
p
p(θ) v det(I(θ)).
Cette distribution est utile dans le cas des variables générées suivant une loi de famille exponen-
tielle, telle que la loi normale. Cette famille de loi se caractérise par une densité de probabilité de la
forme suivante :
avec µ(θ) est une distribution uniforme, c constante de normalisation et m, s sont deux hyperpara-
mètres à déterminer.
Finalement, par l’application du théorème de Bayes, on trouve la loi à posteriori, tel que :
89
9.1.3 Détermination de la distribution à posteriori
Généralement, le point de départ de l’approche bayesienne consiste à définir α(δ) la ddp à priori du
paramètre inconnu f (y1 , , yn /δ) la ddp des observations concernant la variable endogène conditionnée
par δ.
La théorème de bayes permet de calculer la distribution à posteriori du paramètre δ :
f (δ, y1 , yn )
α(δ|y1 , yn ) = ou β(y1 , yn )) est la densité marginale
β(y1 , yn ))
Par ailleurs f (δ, y1 , yn )) = f (y1 , yn ), δ)α(δ)
Ce qui donne α(δ/y1 , yn )) = f (y1 , yn )/δ)α(δ)
Une fois que la distribution à posteriori est déterminée, l’estimateur de δ à posteriori est souvent pris
comme étant égal à l’espérance mathématique de delta ˆ = E(α(δ/y1 , yn )
90
mais qui peuvent ne pas être directement estimées comme faisant partie du modèle original.
Étant donné les avantages d’une approche bayésienne de l’inférence, la question clé est alors la
suivante : à quel point est-il difficile d’intégrer une distribution a posteriori pour produire des résumés
de paramètres ?
Pour certaines distributions, les intégrales pour résumer les distributions postérieures ont des
solutions de forme fermée et sont connues, ou peuvent être facilement calculées en utilisant des
méthodes numériques. Cependant, pour de nombreuses distributions, en particulier multivariées, les
intégrales ne peuvent pas être facilement calculer. Par exemple, si nous avions une distribution a
priori bêta sur la variance d’une distribution normale, la distribution a posteriori de la variance
n’aurait pas de forme connue. Afin de remédier à ce problème, les bayésiens travaillent souvent avec
des a priori conjugués, comme nous l’avons vu dans le chapitre précédent. Cependant, parfois, les
a priori conjugués sont irréalistes, ou un modèle peut impliquer des distributions qui ne sont tout
simplement pas accessibles au calcul simple des quantiles et autres quantités.Dans ces cas, il existe
essentiellement une approche de base pour le calcul des intégrales : les méthodes d’échantillonnage.
La logique de l’échantillonnage est que nous pouvons générer (simuler) un échantillon de taille n à
partir de la distribution d’intérêt et ensuite utiliser des formules discrètes appliquées à ces échantillons
pour approximer les intégrales d’intérêt. Dans une approche d’échantillonnage, nous pouvons estimer
une moyenne en : Z
1X
xf (x)dx = x
n
Différents quantiles peuvent être calculés empiriquement en notant la valeur de x pour laquelle Q%
des valeurs échantillonnées tombent en dessous de celui-ci. Ainsi, l’inférence bayésienne moderne im-
plique typiquement :
(1) établir un modèle et obtenir une distribution a posteriori pour le ou les paramètres d’intérêt,
(2) générer des échantillons à partir de la distribution postérieure et
(3) utiliser des formules discrètes appliquées aux échantillons du distribution a posteriori pour résu-
mer notre connaissance des paramètres.
Ces résumés ne sont pas limités à une seule quantité mais sont virtuellement illimités. Toute sta-
tistique récapitulative que nous calculons couramment pour décrire un échantillon de données peut
également être calculée pour un échantillon à partir d’une distribution a posteriori et peut ensuite
être utilisée pour la décrire.
91
deux méthodes d’échantillonnage, dont chacune est importante pour une compréhension de base des
méthodes MCMC.
où L est la limite inférieure de la densité f. Autrement dit,u = F (z). Donc, exprimé en termes de z :
z = F −1 (u).
Pour donner un exemple concret, prenez la fonction de densité linéaire
1
f (x) = (2x + 3) avec 0 < x < 5
40
il n’y a pas de routines qui permettent d’échantillonner à partir de cette densité, et donc, si l’on
avait besoin de cette densité, il faudrait en développer une. Afin de générer un tirage à partir de
cette distribution en utilisant la méthode d’inversion, nous devons d’abord dessiner u ∼ U (0, 1), puis
calculer z qui satisfait Z z
1
u= (2x + 3)dx
0 40
Nous pouvons résoudre cette équation pour z comme suit. D’abord, évaluez l’intégrale :
92
côté droit de la figure montre les fonctions de densité simulées et théoriques. Remarquez comment
les échantillons des deux densités suivent étroitement, mais ne correspondent pas exactement, les
densités théoriques. Cette erreur est une erreur d’échantillonnage qui diminue à mesure que la taille
de l’échantillon de simulation augmente.
Le programme R suivant a été utilisé pour générer ces tirages.
u=runif(1000,min=0,max=1)
z=(1/2) * (-3 + sqrt(160*u +9))
La première ligne simule 1000 tirages aléatoires de la distribution U (0, 1). La deuxième ligne génère
le vecteur z comme l’inverse de u.
Bien que la méthode d’inversion soit très efficace et facile à mettre en œuvre, deux limitations clés
réduisent sa facilité d’utilisation en tant que méthode générale pour dessiner des échantillons à partir
de densités postérieures. Premièrement, si la fonction inverse est impossible à dériver analytiquement,
la méthode ne peut évidemment pas être utilisée. Par exemple, l’intégrale normale ne peut pas être
résolue directement et, par conséquent, la méthode d’inversion ne peut pas être utilisée pour simuler
à partir de la distribution normale. Dans une certaine mesure, ce problème soulève la question
suivante : , alors pourquoi s’embêter avec la simulation ? Cette question sera abordée sous peu, mais
la réponse courte est que nous ne sommes peut-être pas capables d’effectuer une intégration sur une
densité multivariée, mais nous pouvons souvent décomposer une densité multivariée en univariée pour
laquelle l’inversion peut fonctionner.
Le deuxième problème avec la méthode d’inversion est que la méthode ne fonctionnera pas avec
les distributions multivariées, car l’inverse n’est généralement pas unique au-delà d’une dimension.
Par exemple, considérons la fonction de densité planaire bivariée :
1
f (x, y) = (2x + 3y + 2)
28
avec 0 < x, y < 2. Si on tire u ∼ U (0, 1), nous attendons résoudre double intégrale pour x et y , nous
obtenons :
3
28u = yx2 + xy 2 + 2xy
2
qui, bien sûr, a une infinité de solutions (une équation avec deux inconnues). En pensant à l’avance,
nous pourrions sélectionner une valeur pour une variable, puis utiliser la méthode d’inversion pour
tirer de la distribution conditionnelle de l’autre variable. Ce processus permettrait de réduire le
problème à l’échantillonnage à partir de distributions conditionnelles univariées, ce qui est l’idée de
base de l’échantillonnage de Gibbs, comme je le discuterai bientôt.
93
f (z)
2. Calculez le rapport R = .
m ∗ g(z)
3. Échantillonnez u ∼ U (0, 1). Si R > u, alors accepter z comme un tirage de f (x). Sinon, revenez à
l’étape 1.
Dans cet algorithme, m ∗ g(x) est appelé une "fonction d’enveloppe", en raison de l’exigence que
la fonction de densité g (x) multipliée par une constante m soit supérieure à la valeur de densité pour
la distribution d’intérêt [f ( x)] au même point pour tous les points. En d’autres termes, m ∗ g(x)
enveloppe f (x). À l’étape 1, nous échantillonnons un point z du pdf g (x). Dans l’étape 2, nous
calculons le rapport de la fonction d’enveloppe [m ∗ g(x)] évaluée en z à la fonction de densité
d’intérêt [f (x)] évaluée au même point. Enfin, à l’étape 3, nous tirons U (0, 1) une variable aléatoire u
et la comparons avec R. Si R > u, alors nous traitons le tirage comme un tirage de f (x). Sinon, nous
rejetons z comme venant de f (x), et nous répétons le processus jusqu’à obtenir un tirage satisfaisant.
Cette routine est facile à mettre en œuvre, mais il n’est pas immédiatement évident pour que cela
fonctionne. Examinons à nouveau la densité discutée dans la section précédente et considérons une
fonction d’enveloppe qui est une densité uniforme sur l’intervalle [0, 5] multiplié par une constante
de 2. Je choisis cette constante parce que la hauteur de la densité U (0, 5) est 0.2, alors que la hauteur
maximale de la densité f (x) = (1/40)(2x + 3) est 0.325. Multiplier la densité U (0, 5) par deux
augmente la hauteur de cette densité à 0.4, ce qui est bien au-dessus du maximum pour f (x) et fait
donc de m ∗ g(x) une véritable fonction d’enveloppe. La figure 3 montre les fonctions de densité et
d’enveloppe et représente graphiquement le processus d’échantillonnage de rejet.
Dans la première étape, lorsque nous échantillonnons à partir de la fonction d’enveloppe, nous
choisissons un emplacement sur l’axe des x dans le graphique (voir le graphique du haut dans la
figure 3). Le processus de construction du rapport R et de comparaison avec un écart uniforme est
essentiellement un processus de localisation d’un point dans la direction y une fois la coordonnée
x choisie, puis de décider si elle est sous la densité d’intérêt. Cela devient plus apparent si nous
réarrangeons le rapport et l’inégalité avec u :
f (z)m × g(z) × u
m×g(z)×u nous donne un point dans la dimension y qui se situe quelque part entre 0 et m×g(z).
Ceci peut être facilement vu en notant que m × g(z) × u fournit simplement un tirage au sort à partir
de la distribution U (0, g (z)) : La valeur de ce calcul quand u = 0 est 0 ; sa valeur lorsque u = 1
est m × g(z) (voir le graphique du milieu à la figure 3). Dans la dernière étape, dans laquelle nous
décidons d’accepter z comme un tirage de f (x), nous déterminons simplement si la coordonnée y
tombe en dessous de la courbe f (x) (voir le graphique du bas dans la figure 3). Une autre façon de
penser à ce processus est que le rapport nous indique la proportion de fois où nous accepterons un
tirage à une valeur donnée de x comme provenant de la densité d’intérêt.
Le programme R suivant simule 1000 tirages à partir de la densité f (x) = (1/40)(2x + 3) en
utilisant l’échantillonnage de rejet. La routine tient également compte du nombre de tirages totaux
à partir de g (x) pour obtenir 1000 tirages de f (x).
count=0 ; k=1 ; f=matrix(NA,1000)
94
while(k<1001)
{
z=runif(1,min=0,max=5)
r=((1/40)*(2*z+3))/(2*.2)
if(r>runif(1,min=0,max=1))
{f[k]=z ; k=k+1}
count=count+1
}
La figure 4 présente les résultats d’une exécution de cet algorithme. L’histogramme de l’échan-
tillon de 1 000 correspond très étroitement à la densité d’intérêt. L’échantillonnage par rejet est une
méthode puissante d’échantillonnage à partir de densités pour lesquelles l’échantillonnage par inver-
sion ne fonctionne pas. Il peut être utilisé pour échantillonner à partir de n’importe quelle densité,
et il peut être utilisé pour échantillonner à partir de densités multivariées. Dans le cas multivarié,
nous choisissons d’abord un X - maintenant un vecteur aléatoire, plutôt qu’un seul point - à partir
d’une fonction enveloppante multivariée, puis nous procédons comme précédemment.
L’échantillonnage par rejet a certaines limites. Premièrement, trouver une fonction enveloppante
m × g(z) peut ne pas être une tâche facile. Par exemple, il peut être difficile de trouver une enveloppe
dont les valeurs sont supérieures à tous les points de support pour la densité d’intérêt. Pensez à
essayer d’utiliser une densité uniforme comme une enveloppe pour l’échantillonnage à partir d’une
densité normale. Le domaine de x pour la densité normale va de -1 à +1, mais il n’y a pas de densité
uniforme correspondante. A la limite, une densité U (-1, + 1) aurait une hauteur infiniment faible,
ce qui ferait tomber g (x) en dessous de f (x) au centre de la distribution, indépendamment de la
constante multiple m choisie. Deuxièmement, l’algorithme peut ne pas être très efficace. Si la fonction
enveloppante est considérablement plus élevée que f (x) à tous les points, l’algorithme rejettera la
plupart des tirages tentés, ce qui implique qu’un nombre incroyable de tirages peut être nécessaire
avant de trouver une seule valeur de f (x).
En théorie, l’efficacité d’une routine d’échantillonnage de rejet peut être calculée avant de l’implé-
menter. Dans le cas ci-dessus, la surface totale sous la courbe enveloppante est 2 (5 * 0.4), mais l’aire
totale sous la densité d’intérêt est 1 (par définition d’une fonction de densité). Ainsi, l’algorithme
utilisé devrait accepter environ 50 % des résultats de g (x). En fait, dans le cas présenté et discuté
ci-dessus, il a fallu 2 021 tentatives pour obtenir 1 000 tirages de f (x), ce qui représente un taux de
rejet de 50,5 %. Ces deux limitations rendent l’échantillonnage par rejet, bien que possible, de plus
en plus difficile à mesure que la dimensionnalité augmente dans les distributions multivariées.
95
efficaces, elles ne sont pas encore très utiles en elles-mêmes dans la modélisation statistique du monde
réel. Heureusement, au cours des dernières décennies, des méthodes MCMC ont été développées pour
faciliter l’échantillonnage à partir de distributions complexes.
L’échantillonnage MCMC fournit une méthode pour échantillonner à partir de densités multiva-
riées qui ne sont pas faciles à échantillonner, souvent en décomposant ces densités en des densités
univariées ou multivariées plus gérables. L’approche de base MCMC fournit une prescription pour
(1) échantillonnage à partir d’une ou plusieurs dimensions d’une distribution postérieure et (2) se
déplacer dans tout le support d’une distribution postérieure. En fait, le nom "chaîne de Markov
Monte Carlo" implique ce processus. La partie "Monte Carlo" se réfère au processus de simulation
aléatoire. La partie «Chaîne de Markov» fait référence au processus d’échantillonnage d’une nouvelle
valeur à partir de la distribution postérieure, compte tenu de la valeur précédente : Ce processus
itératif produit une chaîne de valeurs de Markov qui constitue un échantillon de tirages postérieurs.
96
comme fixes.
Nous considérons un exemple simple de Gibbs Sampling basé sur la distribution bivariée :f (x, y) =
1
28
(2x + 3y + 2). La distribution conditionnelle de x est :
f (x, y) 2x + 3y + 2
f (x/y) = =
f (y) 6y + 8
La distribution conditionnelle de y est :
f (x, y) 2x + 3y + 2
f (y/x) = =
f (x) 4x + 10
D ou l’algorithme Gibbs sampling pour échantillonner x et y doit suivre les étapes suivantes :
1. Soit j=0 and fixe les valeurs initiales . Ici soit xj=0 = −5 et y j=0 = −5
2.Tirer xj+1 a partir de f (x/y = y j
3.Tirer y j+1 a partir de f (y/x = xj
4. Incrémenter j = j + 1 et retour à l’étape 2 jusqu’a j = 2000 Comment nous tirons des échantillon à
partir distributions conditionnelles ? Nous savons ce qu’ils sont, mais ce ne sont certainement pas des
distributions standard. Comme ce ne sont pas des distributions standard, mais comme ces conditions
sont univariées et que F −1 () peut être calculé pour chacune d’entre elles, nous pouvons utiliser un
sous-programme d’inversion pour échantillonner chaque densité conditionnelle. Comment trouvons-
nous les inverses de cette densité bivarié ? Rappelez-vous que l’échantillonnage par inversion nécessite
d’abord de dessiner une variable aléatoire U (0, 1), puis d’inverser ce tirage en utilisantF −1 . Ainsi,
pour trouver l’inverse de la densité conditionnelle pour y|x, nous devons résoudre :
Z z
2x + 3y + 2
u=
0 4x + 10
Étant donné que c’est la densité conditionnelle pour y, x est fixe et peut être traitée comme une
constante, et nous obtenons :
Étant donné une valeur courante pour x et un tirage au sort u, z est un tirage au sort de la densité
conditionnelle pour y|x. Un processus similaire peut être entrepris pour trouver l’inverse pour x|y.
Voici un programme R qui implémente l’échantillonnage de Gibbs :
x=matrix(-5,2000) ; y=matrix(-5,2000)
for(i in 2 :2000)
{
# sample from x | y
u=runif(1,min=0, max=1)
x[i]=sqrt(u*(6*y[i-1]+8)+(1.5*y[i-1]+1)*(1.5*y[i-1]+1))-(1.5*y[i-1]+1)
97
#s ample from y | x
u=runif(1,min=0,max=1)
y[i]=sqrt((2*u*(4*x[i]+10))/3 +((2*x[i]+2)/3)*((2*x[i]+2)/3))- ((2*x[i]+2)/3)
}
Ce programme définit d’abord les valeurs de départ pour x et y égales à -5. Ensuite, x est mis à jour
en utilisant la valeur actuelle de y. Ensuite, y est mis à jour en utilisant les valeurs échantillonnées
de x. (Remarquez comment x [i] est calculé en utilisant y [i-1], alors que y [i] est échantillonné en
utilisant x [i].) Les deux sont mis à jour en utilisant la méthode d’échantillonnage d’inversion décrite
ci-dessus. Cet algorithme produit des échantillons à partir des distributions marginales pour x et y,
mais nous pouvons également traiter les paires de x et y comme des tirages de la densité conjointe.
Nous discuterons des conditions dans lesquelles nous pourrons le faire plus en profondeur prochai-
nement. En général, cependant, les distributions marginales pour les paramètres sont particulièrement
intéressantes, car nous nous intéressons souvent à tester des hypothèses concernant un paramètre,
net des autres paramètres d’un modèle. La figure 5 montre un «tracé de traces» des x et des y ainsi
que les densités marginales pour les deux variables. Le tracé de trace est simplement un tracé bidi-
mensionnel dans lequel l’axe x représente l’itération de l’algorithme, et l’axe y représente la valeur
simulée de la variable aléatoire à chaque itération particulière. Heuristiquement, nous pouvons alors
prendre le tracé de la trace, le tourner sur son bord (un tour dans le sens des aiguilles d’une montre)
et laisser l’encre tomber sur l’axe des ordonnées et produire un histogramme du marginal. densité.
Les endroits du tracé qui sont particulièrement foncés représentent des régions de la densité dans
lesquelles l’algorithme simulait fréquemment ; les zones plus légères sont des régions de la densité qui
ont été plus rarement visitées par l’algorithme. Ainsi, "l’encre" s’accumulera plus haut dans les zones
pour lesquelles la variable d’intérêt a une plus grande probabilité. Les histogrammes de ces densités
marginales sont montrés à droite de leurs tracés de trace respectifs, avec les densités marginales
théoriques superposées. Réalisez que ces marginaux ne sont pas normalisés, parce que la constante
de normalisation 1/28 principale annule à la fois le numérateur et le dénominateur.
Notez que, bien que les valeurs initiales aient été très médiocres (-5 n’est pas un point valide dans
les deux dimensions de la densité), l’algorithme a convergé très rapidement vers la région appropriée
[0, 2]. Il faut généralement un certain nombre d’itérations pour un algorithme MCMC pour trouver
la région appropriée et, plus théoriquement, pour la chaîne de Markov produite par l’algorithme
pour échantillonner à partir de la distribution "cible" appropriée. Ainsi, nous rejetons généralement
un certain nombre d’itérations précoces avant de faire des calculs (appelés "burn-in"). Les densités
marginales ne sont donc produites qu’à partir des 1 500 dernières itérations de l’algorithme.
Les histogrammes pour les densités marginales montrent que l’algorithme échantillonne de ma-
nière appropriée à partir des densités d’intérêt. Bien sûr, il y a certainement une erreur - observez
comment les histogrammes ont tendance à être un peu trop bas ou trop haut ici et là. Ceci reflète
une erreur d’échantillonnage, et une telle erreur est réduite en échantillonnant plus de valeurs (par
exemple, en utilisant 5 000 tirages au lieu de 2 000) ;
En plus d’examiner les distributions marginales pour x et y, nous pouvons également examiner la
densité conjointe. La figure 4.6 montre un tracé de trace bidimensionnel, pris à plusieurs étapes. La
98
figure en haut à gauche montre l’état de l’algorithme après 5 itérations ; la figure en haut à droite
montre l’état après 25 itérations ; La figure gauche la montre après 100 itérations ; et la figure en
bas à droite le montre après les 2 000 itérations. Ici encore, nous voyons que l’algorithme, bien que
commençant par de mauvaises valeurs de départ, converge rapidement vers la région planaire partielle
et bidimensionnelle appropriée représentée par f (x, y).
Après échantillonnage de la distribution pour x et y, nous pouvons maintenant résumer notre
connaissance de la densité. La moyenne théorique pour x peut être trouvée en prenant le marginal
pour x (f (x) = (1/28)(4x + 10)) et en intégrant toutes les valeurs pour x :
Z 2
mx = x ∗ f (x)dx = 1.095
0
Un calcul similaire pour y donne une moyenne théorique de 1,143. Les estimations empiriques des
moyennes, basées sur les 1500 derniers tirages des distributions marginales pour les variables (en
éliminant les 500 premières en tant que burn-in) sont x̄ = 1.076 et ȳ = 1.158 . L’écart entre les moyens
théoriques et empiriques est attribuable à l’erreur d’échantillonnage dans l’algorithme MCMC. Une
course plus longue réduirait l’erreur, bien que, même avec 1 500 tirages simulés, les écarts soient
minimes (moins de 2 % pour les x et les y).
99
rapidement vers le support principal du postérieur, conduisant à des temps d’exécution extrêmement
longs.
À l’étape 2, une valeur candidate pour le paramètre ( δ c ) est obtenue en simulant une valeur à
partir d’une densité de proposition [α(.)]. La valeur simulée est considérée comme un «candidat»,
car elle n’est pas automatiquement acceptée comme un tirage de la distribution d’intérêt ; il doit être
évalué pour l’acceptation tout comme dans l’échantillonnage de rejet. En effet, cette étape est un peu
comme si l’on tirait un candidat d’une fonction enveloppante dans l’échantillonnage de rejet, sauf
que la densité de proposition n’a pas besoin d’envelopper réellement la distribution postérieure. Au
lieu de cela, une densité de proposition peut prendre n’importe quelle forme facile à échantillonner
(normale, uniforme, etc.). Étant donné que la densité de la proposition n’est pas la densité d’intérêt,
nous devons vérifier si le paramètre candidat peut être considéré comme provenant de la distribution
cible, tout comme nous l’avons fait dans l’échantillonnage de rejet.
Souvent, nous pouvons utiliser une densité de proposition symétrique (par exemple, une densité
normale ou uniforme) centrée sur la valeur courante du paramètre (δ j−1 ). Par exemple, en utilisant
une densité de proposition normale, le candidat serait tiré d’une distribution normale avec une
moyenne égale à δ j−1 et une certaine variance : δ c = δ j−1 + N (0, c), où c est une constante (plus sur
le choix de c plus tard). Cette approche produit un algorithme de "marche aléatoire de Metropolis",
qui est l’algorithme principal en raison de sa simplicité et de son efficacité. Exemple Estimation
des paramètres normaux via Metropolis.
Considérons n = 1000 valeurs xi générées aléatoirement à partir d’une distribution N (3 ; 25), une
normale avec moyenne m = 3 et variance σ 2 = 25. En utilisant le x généré, nous cherchons à estimer
la moyenne et la variance, les traitant maintenant comme des inconnues. Le plan d’échantillonnage
utilisé dans le code R ci-dessous réellement mises à jour σ plutôt que σ 2 : Définir le vecteur des
inconnues comme δ = (µ, σ 2 ) ; la probabilité est
n
(xi − δ)2
Y 1
p(y/δ) = √ exp −
i=1
( 2πσe ) 2σ 2
1
Supposons un appartement antérieur pour µ (uniforme entre -1 et +1) ; et un avant p(σ) = C’est
σ
a standard noninformative prior on , en fait applée the Jeffreys prior. Ensuite, la densité postérieure
complète est proportionnelle à
100
pratiquement tous les modèles paramétriques en sciences sociales aujourd’hui. Dans cette section, je
me concentre principalement sur la démonstration de plusieurs approches MCMC qui peuvent être
utilisées pour estimer les paramètres du modèle. Le modèle de régression OLS est généralement repré-
senté de deux façons, l’une impliquant une spécification directe d’une distribution pour la variable de
résultat y et l’autre impliquant une spécification de la distribution de l’erreur e. Sous la spécification
classique et typique des sciences sociales, nous supposons que yi est égal à une combinaison linéaire
d’un ensemble de prédicteurs, XiT plus l’erreur ei , et que le terme d’erreur est normalement distribué
avec une moyenne de 0 et une certaine variance σe2 . Sous forme matricielle pour un échantillon entier :
Y = Xβ + e e ∼ N (0, σe2 In )
où In est une matrice d’identité n-dimensionnelle. Souvent, ce terme est omis dans la spécification,
mais la distribution pour le vecteur est techniquement multivariée normale avec un vecteur moyen
0 et une matrice de covariance n × n. Les éléments diagonaux de cette matrice sont tous égaux
(représentant l’homoscédasticité de l’hypothèse des erreurs), et les éléments hors diagonale de cette
matrice sont 0 (représentant l’indépendance de l’hypothèse des erreurs). L’hypothèse de normalité
pour le terme d’erreur n’est pas nécessaire pour l’estimation MCO, mais elle est nécessaire pour
l’estimation du maximum de vraisemblance et pour les tests statistiques classiques sur les paramètres
(les tests t). Sous l’approche du maximum de vraisemblance, une hypothèse d’erreur normale donne
la fonction de vraisemblance suivante :
n
Y
2 −1/2 −1
L(β, σe2 /X, Y
)= (2πσe ) exp T
(yi − Xi β)2
i=1
2σe2
2 −n/2 −1 T
= (2πσe ) exp (Y − Xβ) (Y − Xβ)
2σe2
La solution du maximum de vraisemblance pour trouver les meilleures estimations de β et σe2 (β̂
et σ̂e2 , respectivement) et leurs erreurs-types peut être trouvée en prenant les dérivées première et
seconde du log de cette fonction de vraisemblance. Après dérivation, nous trouvons :
β̂ = (X T X)−1 (X T Y )
1 T
σ̂e2 = e e
n
où e est le vecteur d’erreurs obtenu sous β̂ . Les erreurs types estimées pour les paramètres
peuvent être trouvées en enracinant les éléments diagonaux la matrice de covariance asymptotique
des paramètres, et l’erreur-type pour la variance d’erreur peut également être trouvée :
Plutôt que de spécifier une distribution normale sur le terme d’erreur, une spécification bayésienne
commence généralement par une hypothèse de normalité sur y|x (souvent avec le conditionnement
supprimé) :yi ∼ N (XiT β, σe2 ). Cette spécification donne la même fonction de vraisemblance que la
solution classique.
101
Ce qui reste à faire du modèle entièrement bayésien est la spécification d’un a priori pour β et
σe2 ce qui est souvent fait en spécifiant des a priori indépendants pour chaque paramètre. Un pré-
uniforme incorrect sur la ligne réelle est souvent spécifié pour les paramètres de régression, tandis que
la référence commune avant pour le modèle de distribution normale 1/σe2 est souvent spécifiée pour
le paramètre de variance d’erreur. Cela donne une distribution postérieure qui apparaît comme :
2 2 −n/2+1 −1 T
P (β, σe /X, Y ) = (2πσe ) exp (Y − Xβ) (Y − Xβ)
2σe2
Après simplification. Notez que ce postérieur ne diffère de la fonction de vraisemblance que dans
l’exposant principal. La valeur absolue de l’exposant pour σe2 est augmentée de n / 2 à n / 2 +
1, ce qui est une modification asymptotiquement non pertinente de la fonction de vraisemblance.
Ce résultat souligne une fois de plus que, avec de grands échantillons, le préréglage peut avoir très
peu d’incidence sur l’inférence postérieure. Étant donné une distribution a posteriori, l’objectif de
l’approche bayésienne est de produire plus qu’une simple estimation ponctuelle du paramètre et de
son erreur-type. Ainsi, je discute maintenant plusieurs stratégies pour l’échantillonnage de cette dis-
tribution postérieure. Bien que ce modèle ait été bien étudié, et que les conditions complètes soient
bien connues pour les paramètres, je vais développer un algorithme MH avec deux échantillonneurs
Gibbs différents.
ˆ
richesse i = 0.873 revenui + 0.095agei + −5.795taillei
(34.589) (2.871) (15.597)
102
Intercepter la richesse β̂0 (milliers de dollars), en moyenne, lorsque le revenu, l’âge et la taille de la
famille sont égaux à zéro , on élimine l’intercepte parce qu’elle n’a pas de sens
Coefficients de pente :
β̂1 = 0.873 lorsque le revenu augmente d’une unité (mille dollars), la richesse augmente, en moyenne,
de 0,873 unité (973 dollars) ; en supposant que l’âge et la taille de la famille n’a pas changé.
β̂2 = 0.095 lorsque l’age augmente d’une unité (un ans), la richesse augmente, en moyenne, de 0,095
unité (95 dollars) ; en supposant que l’âge et la taille de la famille n’a pas changé.
β̂3 = −5.795 lorsque le nombre de personnes vivant dans la famille augmente d’une unité (personne),
la richesse diminue en moyenne de 5.795 unités (5795 dollars) ; tenir d’autres facteurs constants.
R2 = 0.2126 le modèle explique 21.26 % de la variation de la richesse financière nette totale
b- Estimation d’un modèle de régression linéaire par l’algorithme Metropolis.
103
# Establish x and y
x=as.matrix(cbind(revenu,age,taille))
y=as.matrix( richesse)
m=2000 # number of iterations
# Establish parameter vectors, proposal scales and acceptance rates s2=matrix(1,m) ;
b=matrix(0,m,3)
bscale=sqrt(diag(vcov(lm(y∼x-1))))*.5
s2scale=sqrt(var(residuals(lm(y∼x-1))*(9275-1)/(9275-3)))*.5
accrate=matrix(0,m,3) ; s2accrate=matrix(0,m)
# unnormalized posterior distribution function
post<-function(x,y,b,s2)
{return((-4636.5*log(s2)+(-.5/s2)*(t(y- x%*%b)%*%(y- x%*%b))))}
# Begin MH Sampling
for(i in 2 :m){
# temporarily set ’new’ values of b
b[i,]=b[i-1,]
# update regression parameters
for(j in 1 :3){
# generate candidate and assume it will be accepted...
b[i,j]=b[i-1,j]+rnorm(1,mean=0, sd=bscale[j]) ; acc=1
#...until it is evaluated for rejection
if((post(x,y,b[i,],s2[i-1]) - post(x,y,b[i-1,],s2[i-1])) <log(runif(1,min=0,max=1)))
{b[i,j]=b[i-1,j] ; acc=0}
accrate[i,j]=(accrate[i-1,j]*(i-1)+acc)/i }
# update s2. generate candidate and assume accepted
s2[i]=s2[i-1]+rnorm(1,mean=0, sd=s2scale) ; acc=1
# ...until it is evaluated for rejection
if(s2[i]<0
(post(x,y,b[i,],s2[i]) - post(x,y,b[i,],s2[i-1])) <log(runif(1,min=0,max=1)))
s2[i]=s2[i-1] ; acc=0 s2accrate[i]=(s2accrate[i-1]*(i-1)+acc)/i
if(i%%10==0){print(c(i,b[i,1],b[i,2],b[i,3],s2[i],accrate[i,1],s2accrate[i]))} }
summary(b[,1])
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8535 0.8732 0.8334 0.8923 0.9465 summary(b[,2])
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.008471 0.063852 0.085447 0.0975 0.1143 0.364149 summary(b[,3])
Min. 1st Qu. Median Mean 3rd Qu. Max.
-6.720 -5.957 -5.682 -5.395 -5.383 1.202
104
9.5.2 Gibbs sampling
Il y a au moins deux façons générales de développer un échantillonneur de Gibbs pour le modèle
de régression linéaire. La première méthode consiste à déterminer les conditions complètes pour (1) le
vecteur du paramètre de régression et (2) le paramètre de variance d’erreur. La distribution complète
postérieure conditionnelle pour le paramètre de variance d’erreur est simple à dériver de l’équation
9. Avec β connu / fixé, le postérieur conditionnel pour σe2 est :
T
2 2 −n/2+1 −e e
P (σe /β, X, Y ) = (σe ) exp )
2σe2
où eT e est la somme des termes d’erreur carrés sous la valeur donnée pour β . Ce postérieur condi-
tionnel est facilement considéré comme une distribution gamma inverse avec les paramètres α = n2 et
eT e
β= . La distribution conditionnelle pour β est légèrement plus difficile à dériver, car sa dérivation
2
implique une algèbre matricielle. Cependant, le processus est identique dans le concept à l’approche
que nous avons utilisée dans l’exemple de distribution normale univariée . Avec σe2 fixe, nous pouvons
nous concentrer exclusivement sur l’exponentielle :
−1 T
exp (Y − Xβ) (Y − Xβ)
2σe2
Premièrement, nous pouvons distribuer la transposition dans le premier terme du numérateur :
−1 T T T
exp (Y − β X )(Y − Xβ)
2σe2
puis développons la multiplication :
−1 T T T T T T
exp (Y Y − Y X − β X Y + β X Xβ)
2σe2
Le premier terme est constant par rapport à β, et peut donc être éliminé en tant que constante de
proportionnalité multiplicative. Les deux termes du milieu sont identiques l’un à l’autre - l’un est
juste une version transposée de l’autre, mais les deux sont 1.1 ; ainsi, on peut être transposé, et les
deux peuvent être groupés. Après avoir réarrangé les termes, nous obtenons :
−1 T T T T
exp (β X Xβ + 2β X Y )
2σe2
105
sa distribution normale multivariée avec σe2 fixe, et (3) échantillonnant σe2 de son distribution inverse
gamma avec β (et donc e) fixe.
Estimation d’un modèle de régression linéaire par l’algorithme Gibbs sampling.
m=5000 # number of iterations
# establish parameter vectors and constant quantities
s2=matrix(1,m) ; b=matrix(0,m,3)
xtxi=solve(t(x)%*%x)
pars=coefficients(lm(y∼x-1))
# Gibbs sampling begins
for(i in 2 :m){
# simulate beta from its multivariate normal conditional
b[i,]=pars+t(rnorm(3,mean=0,sd=1))%*%chol(s2[i-1]*xtxi)
# simulate sigma from its inverse gamma distribution
s2[i]=1/rgamma(1,9275/2,.5*t(y-x%*%(b[i,]))%*%(y-x%*%(b[i,])))
# write output to file and screen
if(i%%50==0)print(c(i,b[i,1],b[i,2],b[i,3],s2[i]))
}
summary(b[,1])
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8576 0.8735 0.8738 0.8904 0.9588
summary(b[,2])
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.03135 0.07194 0.09440 0.09438 0.11698 0.21624
summary(b[,3])
Min. 1st Qu. Median Mean 3rd Qu. Max.
-7.017 -6.046 -5.796 -5.794 -5.542 0.000
106
Chapitre 10
107
Xt−1 = LXt
Xt−2 = L2 Xt = L(LXt )
.
.
Xt−k = Lk Xt
De même un polynôme de retard est un polynôme en L d’ordre fini ou infini :
φ(L) = φ0 + φ1 L + ........ + φp Lp = pi=1 (φi Li ) avec « p = ordre fini »
P
Yt = ∆Xt = Xt − Xt−1
set.seed(123)
x <- as.ts(rnorm(6)) # génération de 6 valeurs aléatoires normales de classe ts.
lag1.x <- lag(x, -1) # série retardée d’une seule observation.
lag2.x <- lag(x, -2) # série retardée de deux observations.
result <- t(cbind(x, lag1.x, lag2.x))
res2 <- ts.intersect(x, lag1.x, lag2.x)
deltaX <- diff(x) # première différence
delta2X <- diff(xt, differences = 2) # différence d’ordre 2
108
cov(xt − x∗t , xt−k − x∗t−k ) cov(xt − x∗t , xt − x∗t−k )
rp (k) = p =
ν(xt − x∗t )ν(xt−k − x∗t−k ) ν(xt − x∗t )
Où x∗t (respectivement x∗t−k ) est la régression affirme de xt (respectivement x∗t−k ) sur xt−1 .......xt−k−1
Cette expression montre que rp (k) est le coefficient de régression de : (xt − x∗t )sur (xt−k − x∗t−k )
En langage R : Les fonctions acf et la pacf calculent la corrélation et la corrélation partiellle.
cov(xt , xs ) τ (t, s)
ρ(t, s) = p =p
cov(xt , xt )cov(xs , xs ) τ (t, t)τ (s, s)
109
Dans le cas d’un processus stationnaire on calculera simplement :
τ (t − s) τ (t)
ρ(t, s) = = = ρ(τ )oτ = t − s
τ (0) τ (0)
Lorsque le processus stochastique est stationnaire, on peut dessiner le graphique des différents co-
efficients d’autocorrélation en fonction du retard (ou écart de dates) h. Ce graphique est l’auto-
corrélogramme théorique du processus. Il indique comment une composante de date quelconque du
processus est liée linéairement aux composantes des dates précédentes, pour des retards croissants.
C’est une représentation graphique de la mémoire du processus, qui montre dans quelle mesure ses
réalisations courantes sont influencées par ses réalisations passées.
E(yt ) = E(yt−s ) = µ
V (yt ) = V (yt−s ) = σ 2
Cov(yt , yt−s ) = Cov(yt−j , yt−j−s ) = γ
où µ, σ 2 et γ sont constantes
Ainsi une série temporelle est dite non stationnaire si au moins un de ces deux premiers moments
évolue avec le temps.
Exemple d’une série stationnaire
set.seed(123)
bb = rnorm(500, 0, 1)
plot.ts(bb, xlab = "", ylab = "", panel.first = grid(col = "blue"), col = "red")
Ce processus peut être exprimer comme une fonction déterministe du temps plus un processus
stationnaire d’espérance mathématique nulle et de variance constante. Une série non stationnaire
110
déterministe elle s’écrit sous la forme suivante :
yt = µ + βt + ut
E(yt ) = µ + βt
V (yt ) = σ 2
L’espérance de yt dépend donc du temps, mais sa variance ne dépend pas du temps. Ce type de
processus est stationnaire uniquement en écart par rapport à une tendance déterministe et dite Ts
(Trend stationnaire)
Ce processus est caractérisé par la présence d’au moins une racine unitaire. Une série dont la non
stationnarité, est due à la présence d’une racine unitaire, elle s’écrit sous la forme suivante :
yt = µ + yt−1 + νt
Quel on peut écrire sous la forme suivante :
t−1
X
yt = µt + νt−1 + y0
i=1
E(yt ) = y0 + µt
V (yt ) = tσ 2
On voit ici que la moyenne ainsi que la variance sont en fonction du temps. Il est stationnaire
uniquement qu’après la différenciation de la série.
Exemples des processus non stationnaires
X1 <- cumsum(rnorm(500)) # marche aléatoire.
X2 <- c(rnorm(200, 0, 10), 40 + rnorm(300, 0, 10)) # processus avec rupture.
plot.ts(X1, xlab = "", ylab = "", main = "Processus marche aléatoire")
plot.ts(X2, xlab = "", ylab = "", main = "Processus avec rupture")
111
de cette série.
On souhaite estimer un modèle linéaire reliant plusieurs séries économiques (par exemple une
variable dépendante et une ou plusieurs variables explicatives) et tester des hypothèses sur les para-
mètres de cette relation. Les techniques d’inférence statistique à utiliser sont différentes en fonction
de la nature des processus stochastiques de ces séries :
Les techniques d’inférence statistique classique sont valables à condition qu’aucune série de la re-
lation n’ait une tendance stochastique. Les techniques liées à la problématique de la cointégration
s’imposent lorsque des séries de la relation ont une tendance stochastique.
Les tests de racine unitaire permettent de tester l’hypothèse H0 qu’une série économique observée
est la réalisation d’un processus stochastique non stationnaire à tendance stochastique, contre l’hy-
pothèse H1 que ce processus est stationnaire à tendance uniquement déterministe.
Ces tests sont basés sur l’estimation préalable d’un modèle autorégressif AR(p) censé approcher le
vrai comportement du processus stochastique dont la série observée est une réalisation. Un test de
non stationnarité largement utilisé et répandu est le test de racine unitaire proposé par Dickey et
Fuller en 1979. L’hypothèse nulle du test est la présence de racine unitaire, soit le non stationnarité
de type stochastique. Le test consiste à tester :
H0 : φ = 1 contre H1 : φ < 1
dans le modèle :
Yt = φYt−1 + εt
avec εt bruit blanc (0, σ 2 ). L’hypothèse nulle correspond au cas de marche aléatoire pure (processus
DS, I(1)) et l’hypothèse alternative correspond au cas d’un modèle AR(1) stationnaire.
Pour mener ce test, on calcule la statistique de Student, mais attention, cette statistique ne suit
plus sous l’hypothèse nulle une loi de Student, puisque, sous l’hypothèse nulle, le processus est non
stationnaire de type DS et les propriétés asymptotiques ne sont plus standards. Ainsi, la différence
avec un test standard repose sur les valeurs critiques à utiliser pour conclure sur le test. On ne
peut plus utiliser 1.64 comme valeur critique pour un test unilatéral à 5%. Il faut utiliser les valeurs
critiques, qui ont été retabulèes par Dickey et Fuller.
De plus, ce test ne répond pas a nos attentes de détection du type de non stationnarité dans les
variables économiques, d’une part parce que l’hypothèse de processus TS n’est pas présente et d’autre
part parce que les séries économiques sont caractérisées par de l’auto corrélation, qui conduira la
plupart du temps à rejeter l’hypothèse de bruit blanc pour εt dans le test ci-dessus. Pour prendre
en compte, d’une part la présence d’auto corrélation dans les sériés économiques, et, d’autre part,
l’hypothèse de tendance déterministe, on mène les tests de racine unitaire dans les trois régressions
suivantes. p
X
∆Yt = ρYt−1 + α + β.t + φj ∆Yt−j + εt (1)
j=1
p
X
∆Yt = ρYt−1 + α + φj ∆Yt−j + εt (2)
j=1
112
p
X
∆Yt = ρYt−1 + φj ∆Yt−j + εt (3)
j=1
avec p le nombre de retards à ajouter dans la régression afin de prendre en compte l’auto corré-
lation et donc de “blanchir” les résidus. On parle de correction paramétrique de l’auto corrélation et
on appelle les tests de Dickey-Fuller, les tests de Dickey-Fuller augmenté (ADF).
H0 : ρ = 1 contre H1 : ρ < 1
• Dans le modèle (1) : Yt est I(1) + T sous H0 .
Sous H1 , Yt a une tendance déterministe et l’écart à cette tendance déterministe suit un modèle AR
stationnaire, on note I(0) + T, soit un processus TS.
• Dans le modèle (2) : Yt ,sous H0 est I(1) + C .
Sous H1 , Yt suit un modèle AR stationnaire non centré, on note I(0) + C.
• Dans le modèle (3) : Yt sous H0 , est I(1).
Sous H1 , Yt suit un modèle AR stationnaire, on note Yt ∼ I(0).
Les valeurs critiques du test de racine unitaire dépendent de la présence ou non d’une constante ou
d’une tendance. Ainsi, il faut comparer la statistique de student à la valeur critique pertinente. Soit,
au seuil de 5%, on doit considérer comme valeur critique, non pas -1.64 comme dans le cas standard,
mais -3.45 dans le modèle (3), -2.89 dans le modèle (2) et -1.95 dans le modèle (1).
Sous R, les tables sont intégrées dans le logiciel, et R donne directement les p-values calculées selon
les bonnes valeurs critiques. De plus, puisque les valeurs critiques dépendent de la présence ou non
d’une constante ou d’une tendance, cela implique que le test de racine unitaire doit être mené dans
le "bon modèle".
Ainsi, une possibilité de mettre en oeuvre les tests de racine unitaire est de procéder de manière
emboitée, selon la stratégie suivante : on teste la racine unitaire dans le modèle le plus General, puis
on teste si le modèle utilise pour mener le test était pertinent. Si tel n’est pas le cas, on doit mener
à nouveau le test de racine unitaire dans le modèle contraint, etc.
x <- rnorm(100) # no unit-root
y <- diffinv(x) # contains a unit-root
adf.test(x, alternative="stationary", k=0)
adf.test(y, alternative="stationary", k=0)
library("vars")
summary(ur.df(x, type = "trend", lags = 2))
summary(ur.df(y, type = "trend", lags = 2))
113
des composantes courante et passées d’un bruit blanc. La représentation de Wold du processus
stochastique stationnaire Xt se traduit en termes mathématiques de la façon suivante :
Xt = µ + φ1 εt + φ2 εt−1 + ...∀t ∈ Z
εt : un bruit blanc.
φ1 , φ2 , .. des coefficients.
Ici, l’espérance de Xt est représentée par µ et la variance de Xt par σa2 ( pi=1 φ2i + 1)
P
AR(1) : yt = φ1 yt−1 + εt
AR(2) : yt = φ1 yt−1 + φ2 yt−2 + εt
.
.
AR(p) : yt = φ1 yt−1 + ........... + φp yt−p + εt
Où les φi sont des paramètres à estimer et εt est un bruit blanc de variance σ 2 .
Nous pouvons ajouter a ce processus une constante qui ne modifie en rien les propriétés stochas-
tiques.
a-Stationnarité d’un processus autorégressif :
Pour un processus AR (1) on a :
yt = φ1 yt−1 + εt
Ce modèle peut être écrit aussi ainsi :
X
yt = φn1 εt−n + εt−n+1 + . . .
114
Pour un processus AR(p), on a :
p
X
yt = φj yt−j + εt
j=1
Ou encore Φ(L)yt = εt
Avec Φ (L) = 1 − φ1 L − · · · − φp Lp ou L est un opérateur de retard
Φ (L) est appelé le polynôme retard associé à l’AR(p).
Le polynôme caractéristique associé à l’AR(p) est défini par :
Φ (z) =z p − φ1 z p−1 − φ2 z p−2 − · · · − φp−1 z − φp avec L = z1
Les polynômes Φ (L) et Φ (z) sont utilisés pour donner une condition de stationnarité au second ordre
1
du processus AR. Si Le est une racine du polynôme retard alors Lf = z̃ est une racine du polynôme
caractéristique.
Le processus AR(p) est asymptotiquement stationnaire au second ordre si et seulement
si les racines du polynôme retard Φ (L) sont à l’extérieur du cercle unité, c’est-à-dire
plus grande que l’unité en module.
Ou encore,
Le processus AR(p) est asymptotiquement stationnaire au second ordre si et seulement
si les racines du polynôme caractéristique Φ (z) sont à l’intérieur du cercle unité, c’est-
à-dire plus petite que l’unité en module.
Exercice :
On considère le processus AR(2) suivant :
γ0 = φ1 γ1 + φ2 γ2 + · · · + φp γp + σε2
γ1 = cov (yt , yt−1 ) = E (yt yt−1 ) = E ((φ1 yt−1 + φ2 yt−2 + · · · + φp yt−p + εt ) yt−1 )
= φ1 E yt−1 2 + φ2 E (yt−1 yt−2 ) + · · · + φp E (yt−1 yt−p ) + E (yt−1 εt )
γ1 = φ1 γ0 + φ2 γ1 + · · · + φp γp−1
De même, on peut montrer que :
γ2 = φ1 γ1 + φ2 γ0 + · · · + φp γp−2
115
.
γk = φ1 γk−1 + φ2 γk−2 + · · · + φp γp−k
Ces équations sont appelées : les équations de Yule-Walker
On peut déterminer les coefficients φ1 en utilisant ces équations de Yule-Walker. En effet, on prend un
simple exemple pour un AR(2). Les équations de Yule-Walker sont présentées sous la forme suivante :
(
γ1 = φ1 γ0 + φ2 γ1
γ2 = φ1 γ1 + φ2 γ0
La résolution de ce système nous donne
γ1 γ1
γ2 γ0
φ1 =
γ0 γ1
γ1 γ0
γ0 γ1
γ1 γ2
Et φ2 =
γ0 γ1
γ1 γ0
P
t −y)(yt+k −y)
Sachant que la fonction d’autocorrélation ρk = covariancevariance
pour le retard k
= (yP (yt −y)2
= γk
γ0
Les équations de Yule-Walker peuvent être aussi représentées sous cette forme :
ρ1 = φ1 ρ0 + φ2 ρ1 + · · · + φp ρp−1
ρ2 = φ1 ρ1 + φ2 ρ0 + · · · + φp ρp−2
ρk = φ1 γk−1 + φ2 ρk−2 + · · · + φp ρ0
L’autocorrélogramme simple d’un processus AR(p) stationnaire décline rapidement avec l’aug-
mentation du retard, mais il ne s’annule jamais.
La fonction d’autocorrélation partielle d’un processus AR(p) vérifie :
rp (k) = 0 ∀k > p
Un processus AR(p) reste donc lié linéairement à toutes ses valeurs passées, aussi éloignées soient-
elles.
Simulation d’un AR(1) : yt = 0.7yt−1 + εt
ts.ar<- arima.sim(list(order = c(1,0,0), ar =0.7), n = 200)
ts.plot(ts.ar,type="l")
acf(ts.ar,10)
pacf(ts.ar,10)
116
AR1.png
Pour un processus AR(1) stationnaire, pour tout h, et l’autocorrélogramme peut montrer :
-soit une diminution exponentielle sans changement de signe ;
-soit une diminution exponentielle avec changement de signe à chaque retard.
En revanche, l’autocorrélogramme partiel s’annule après le premier retard : r1 (k) = 0∀k > 1
Simulation d’un AR(2) : yt = 0.8yt−1 − 0.5yt−2 + εt
ts.ar<- arima.sim(list(ar = c(0.8, -0.5)), n = 200)
ts.plot(ts.ar,type="l")
acf(ts.ar,10)
pacf(ts.ar,10)
AR2.png
Pour un processus AR(2 ) stationnaire, l’autocorrélogramme simple peut montrer :
-soit une diminution exponentielle sans changement de signe .
-soit une diminution exponentielle avec changement de signe à chaque retard .
-soit une sinusoide s’amoindrissant.
En revanche, l’autocorrélogramme partiel s’annule après le deuxième retard : r2 (k) = 0∀k > 2
117
10.2.3 Processus moyenne mobile d’ordre q
Dans ce processus chaque observation xt est générée par une moyenne pondérée d’aléas jusqu’à
la qème période :
M A(1) : xt = εt − θ1 εt−1
M A(2) : xt = εt − θεt−1 − θ2 εt−2
M A(q) : xt = εt − θ1 εt−1 − ........ − θq εt−q
Où les θi sont les paramètres et εt est un bruit blanc.
L’équation précédente peut s’écrie sous forme :
xt = (1 − θ1 L − θ2 L2 − ...... − θq Lq )εt
Donc, un processus de moyenne mobile n’est qu’une combinaison linéaire de termes d’erreur de nature
bruit blanc. Il est à noter qu’il y a une équivalence entre un processus M A (1) (yt = εt + θ1 εt−1 ) et
un processus AR d’ordre infini :
M A (1) = AR(∞)
Cette équivalence n’est possible que sous l’hypothèse |θ1 | < 1. Cette condition est appelée condition
d’inversibilité du processus.
Ainsi, pour un MA (q) si les racines de Θ(L) sont toutes de module supérieur à 1, on dit
que le processus est inversible.
a- Stationnarité d’un processus à moyenne mobile :
Pour un processus MA(1) on a :
xt = εt − θ1 εt−1
E(y t ) = E (εt − θ1 εt−1 ) = E (εt ) − θ1 E(εt−1 ) = 0
Donc E(y t ) est nulle et égale à zéro.
V ar(y t ) = V ar (εt − θ1 εt−1 ) = V ar (εt ) + θ12 V ar (εt−1 ) = (1 + θ12 )σε2
la variance est aussi constante dans le temps
118
1 + θ12 + θ22 + .. + θq2 σε2 si h = 0
cov(Y t Yt−h ) = (−θh + θ1 θh+1 + · · · + θq−h θq ) σε2 si 0<h<q
0 si h > q
Sont toutes constantes dans le temps, en conséquence, un processus à moyenne mobile est
toujours stationnaire. De même, on peut montrer que la fonction d’autocorrélation d’un processus
MA(q) est égale à :
1 si h = 0
cov(Y t Yt−h ) (−θh +θ1 θh+1 +···+θq−h θq )
ρh = p = si 0<h<q
var(Yt )var(Yt−h ) (1+θ12 +θ22 +..+θq2 )
0 si h > q
Pour un processus MA(1), l’autocorrélogramme simple s’annule après le premier retard (ρh = 0 pour
tout h > 1), tandis que l’autocorrélogramme partiel peut montrer :
-soit une diminution exponentielle sans changement de signe ;
-soit une diminution exponentielle avec changement de signe à chaque retard ;
119
Simulation d’un MA(2) : yt = 0.7εt−1 + 0.4εt−2 + εt
ts.ma2 = arima.sim(list(ma =c(0.7,0.4)), n = 200)
ts.plot(ts.ma,type="l")
acf(ts.ma,10)
pacf(ts.ma,10)
Pour un processus MA (2), l’autocorrélogramme simple s’annule après le deuxième retard (ρh = 0
pour tout h > 2) tandis que l’autocorrélogramme partiel peut montrer :
-soit une diminution exponentielle sans changement de signe
-soit une diminution exponentielle avec changement de signe à chaque retard
-soit une sinusoïde s’amoindrissant
120
Nous avons, ARM A(1, 0) = AR(1) et ARM A(0, 1) = M A(1)
Le processus ARM A(p, q) est stationnaire si et seulement si les racines de polynôme
Φ (L) en module sont toutes supérieur à 1.
Le processus ARM A(p, q) admet une inverse si et seulement si les racines de Θ (L) sont
supérieur à 1.
121
ou encore : Φ (L) yt = Θ (L) ut avec (1 − L)d xt = yt
10.3.2 Identification
La phase d’identification est la plus importante et la plus difficile : elle consiste à déterminer le
modèle adéquat dans la famille des modèle ARIMA. Elle est fondée sur l’étude des corrélogramme
122
simple et partiel. Nous pouvons essayer d’édicter quelques règles simples facilitant la recherche des
paramètres p, d, q du modèle ARIMA.
Après stationnarisation, nous pouvons identifier les valeurs des paramètres p, q du modèle ARM A a
partir de corrélogramme qui est tracé par les fonctions d’auto-corrélation empiriques. Ces fonctions
vérifient les propriétés suivantes :
E(ρ̂k ) = − n1 et varρ̂k ) = n1 où n est la taille de la série en question.
ρ̂k est asymptotiquement normalement distribuées. D’où pour tester :
1. Si la fonction d’autocorrélation simple n’a que ses q premiers termes différents de 0 et que
les termes da la fonction d’autocorrélation partiel diminuent géométriquement, il s’agit d’un
M A(q).
2. Si la fonction d’autocorrélation partiel n’a que ses p premiers termes différents de 0 et que
les termes da la fonction d’autocorrélation simple diminuent géométriquement, il s’agit d’un
AR(p).
3. Si les fonctions d’autocorrélation simple et partiel ne paraissent pas tronquées, il s’agit alors
d’un processus de type ARMA, dont les paramètres dépendent de la forme particulière des
corrélogrammes.
10.3.3 Estimation
Il existe plusieurs méthodes d’estimation des paramètres (MCO, Yule-Walker, MV). La méthode
la plus repandue, est celle du maximum de vraisemblance (MV). Cette méthode repose sur l’hypothèse
de normalité des résidus.( ∼ N (0, σ 2 ) iid . Le logarithme de vraisemblance d’un processus ARMA(p ;
q) est donnée par :
n n 1 0 S(φ, δ)
lnLn = − ln(2π) − ln(σ 2 ) − ln(det(ZZ ) −
2 2 2 2σ 2
Z est une matrice de taille (p + q + n, p + q) qui dépends des paramètres φi , i = 1..p, δj , j = 1..., q.
En maximisant cette fonction, on déduit les estimateurs φ̂i , δˆj et σˆ2
123
10.3.4 Validation
Souvent, il n’est pas facile de déterminer un modèle unique qui représente le processus générateur
de données, et il n’est pas rare pour estimer plusieurs modèles à l’étape initiale. Le modèle qui est
finalement choisi est celui considéré comme le meilleur basé sur un ensemble de critères de contrôle
et de diagnostic. Ces critères comprennent les t-tests de significativité des paramètres estimés.
Les résidus issues d’une estimation doivent vérifier quelques propriétés statistiques :
- La normalité : La normalité peut être tester graphiquement, soit en représentant l’histogramme
des résidus, soit par le graphe quantile-quantile (qq-plot).
- Absence d’auto-corrélation : Il existe plusieurs tests d’absence d’auto-corrélation. Ces tests se re-
groupent en deux groupe. Un groupe regroupe les tests paramétriques et un groupe constitué des
tests non paramétriques. Les tests paramétriques les plus connus dans la littérature étant ceux de
Box et Pierce (1970) et de Ljung et Box (1978).
Le test de Box-Pierce permet d’identifier les processus sans mémoire (suite de variables aléatoires
indépendantes entre elles). Nous devons donc identifier cov(yt , yt−k ) = 0 ou encore ρk = 0∀k.
Un processus de bruit blanc implique que ρ1 = ρ2 = ... = ρh = 0 , soit les hypothèses :
H0 : ρ1 = ρ2 = ... = ρh = 0
H1 : il existe au moins un ρi significativement différent de 0.
Pour effectuer ce test, on recourt à la statistique Q qui est donnée par :
h
X
Q=n ρ2k
k=1
124
On utilise la même statistique que précédemment, c-à-d LM (p) = nR2 Si LM < X1−α/2 2
(2p), on
accepte l’hypothèse nulle
Analyse des résidus
ts.plot(fit4$resid,type="l",col="blue")
abline(0,0,col="red")
acf(fit4$resid[2 :200],10,"correlation") on obtient Box.test(fit4$resid, lag = 10, type =
"Box-Pierce") # absence d’auto-corrélation :
X-squared = 11.9391, df = 10, p-value = 0.2892
Donc les résidus ne sont pas autocorrélés. -Pour le contrôle graphique de la normalité des
résidus :
qqnorm et qqline.
125
Chapitre 11
La modélisation VAR
Dans ce qui a précède, nous avons supposé que l’on s’intéresse à une seule série temporelle. En
pratique, on dispose souvent d’observation simultanées de plusieurs variables : ceci permet d’étudier
les liaisons existant entre elles, de voir comment celle-ci évoluent avec le temps, d’examiner le délai
nécessaire pour qu’une variable influe sur une autre. Il s’agit dans ce suit de généraliser les résultats
présentes dans les cas unidimensionnel.
L’application des vecteurs autorégressifs (VAR) est donc le prolongement de l’analyse dynamique
des phénomènes économiques faisant intervenir nombreuses variables endogènes expliqués par leurs
valeurs retardées et par leurs valeurs retardées d’une autre variable endogène. Les modèles VAR sont
une combinaison des modèles à équations simultanées et des processus autorégressifs, qui permet la
distinction entre variables exogènes et endogènes.
11.1 Définitions
Plusieurs auteurs ont définir un processus VAR comme un ensemble de n variables, yit pour
i=1,2,. . ., n, observées dans le temps t= 1,2,. . ., T chacune de ces variables est liées d’une manière
linéaire, a son propre passé ainsi qu’au passé immédiat et le passé lointain de p périodes de toutes
les autres variables.
126
Pour le cas de deux variables n=2 et un seul retard p=1, la formalisation correspond au système
suivant :
yt = Ayt−1 + ut
Exemple le modèle VAR d’un ordre p = 4
Soit une représentation VAR dans laquelle on considère deux variables y1t et y2t . Chacune de ces
variables est fonction de ses propres valeurs passées et de celles de l’autre. Par exemple le modèle
VAR d’un ordre p = 4, s’écrit :
4
X 4
X
y1t = a1 + b1i y1t−i + c1i y2t−i − d1 y2t + 1t
i=1 i=1
4
X 4
X
y2t = a2 + b2i y1t−i + c2i y2t−i − d2 y1t + 2t
i=1 i=1
Les variables y1t et y2t sont considérées comme étant stationnaires, les perturbations 1t et 2t (les
innovations ou les chocs) sont des bruits blancs de variances constantes σ21 et σ22 et non autocorrélées.
Nous pouvons immédiatement constater l’abondance de paramètres à estimer (ici 20 coefficients) et
les problèmes de perte de degrés de liberté qui en résultent. À la lecture de ce modèle, il apparaît
qu’il n’est pas sous forme réduite : en effet, y1t a un effet immédiat sur y2t et réciproquement y2t a
un effet immédiat sur y1t . Ce système initial est appelé forme structurelle de la représentation VAR.
Sous forme matricielle, ce modèle devient :
4
X
BYt = A0 + Ai y1t−i + t
i=1
" # " # " # " #
1 d1 y1t bi1 c1i 1t
Avec B = Y = ;A = ;=
d2 1 y2t b2i c2i 2t
127
Pour obtenir un modèle VAR sous forme standard, on multiplie de part de d’autre dans (1) par
−1
B . Le modèle sous forme standard s’écrit :
4
X 4
X
y1t = a01 + a11i y1t−i + a21i y2t−i + v1t
i=1 i=1
4
X 4
X
y2t = a02 + a12i y1t−i + a22i y2t−i + v2t
i=1 i=1
Un groupe de variables aléatoires temporelles est généré par un modèle VAR si chacune de ses
variables est une fonction linéaire de ses propres valeurs passées et des valeurs passées des autres
variables du groupe, à laquelle s’ajoute un choc aléatoire de type bruit blanc.
Définition 2 : On applique un modèle VAR dans le cas d’un groupe des variables aléatoires qui
dépend de leur passées ainsi que le passées des autres variables. La généralisation de la représentation
VAR a k variables et p retards s’écrit sous la forme matricielle suivante :
128
les tests de racine unitaire, tels que les tests Dickey-Fuller.
Yt = V + A1 Yt−1 + · · · + Ap Yt−p + εt
Où εt est un bruit blanc de matrice de variance-covariance Ω.
La vraisemblance conditionnellement à toutes les valeurs passées du processus s’écrit :
T
Y
L (yt , . . . , yt−1 ) = L(yt |yt−1 )
t=1
T
Y 1
L (yt , . . . , yt−1 ) = √ N √
t−1 2π detΩ
1X
exp − (Yt − A1 Yt−1 − · · · − Ap Yt−p )Ω−1 (Yt − A1 Yt−1 − · · · − Ap Yt−p )
2
La log-vraisemblance du processus VAR(p), s’écrit donc :
T
TN T 1X
logL (y1 , . . . , yT ) = − log2π − logΩ − εt Ω−1 εt
2 2 2 t=1
La maximisation de cette log-vraisemblance permet alors d’obtenir les estimations des paramètres
A1 , ..., Ap et de la matrice de variance-covariance Ω.
129
Dans le cas d’un processus VAR, chacune des équations peut être estimée par les MCO, indé-
pendamment les unes des autres (ou par une méthode de maximum de vraisemblance). De même,
nous pouvons ajouter à la spécification VAR des variables binaires afin de corriger un mouvement
saisonnier ou une période anormale
T + n∗ K X
F P E(n) = ( ) det ( (n))
T − n∗ u
Avec u (n) =T −1 Tt=1 t 0 t et n∗ est le nombre total de paramètres dans chaque équation et n
P P
attribue l’ordre de décalage.On choisir l’ordre p qui minimise les critères d’informations. Pour dé-
terminer le nombre de retards d’un modèle à retards échelonnés, nous avons présenté les critères de
Akaike et de Schwarz. Dans le cas de la représentation VAR, ces critères peuvent être utilisés pour
déterminer l’ordre p du modèle. La procédure de sélection de l’ordre de la représentation consiste
à estimer tous les modèles VAR pour un ordre allant de 0 à h (h étant le retard maximum admis-
sible par la théorie économique ou par les données disponibles). Les fonctions AIC(p) et SC(p) sont
calculées de la manière suivante :
2k 2 p
AIC(p) = Ln[detΣe ] +
n
k 2 pLn(n)
SC(p) = Ln[detΣe ] +
n
avec : k = nombre de variables du système ; n = nombre d’observations ; p = nombre de retards ; Σe
= matrice des variances covariances des résidus du modèle.
Le retard p qui minimise les critères AIC ou SC est retenu.
130
11.1.3 Application sous R
Nous cherchons à modéliser, sous la forme VAR, la demande (y1t ) et les prix (y2t ) d’une matière
première. Nous disposons des données trimestrielles CVS et en différences premières sur 18 ans .
On demande :
1) de rechercher l’ordre du modèle VAR,
2) d’estimer les paramètres du modèle,
3) de calculer une prévision pour l’année 19 avec son intervalle de confiance à 95 %.
Solution
Nous allons utiliser les critères de Akaike et de Schwarz pour des décalages h allant de 0 à 4. Nous
devons donc estimer quatre modèles différents et retenir celui dont les critères AIC et SC sont les
plus faibles.
library("vars")
tab1=read.table("modeleVAR.txt",h=T)
attach(tab1)
layout(matrix(1 :2, nrow = 1, ncol = 2))
plot.ts( Y1)
plot.ts( Y2)
y=cbind(Y1,Y2)
# On détermine une longueur de retard optimale pour un VAR
VARselect(y, lag.max = 4, type = "both")
131
Y2.l1 0.2993 0.1223 2.446 0.0170 *
const -12.8628 5.3352 -2.411 0.0186 *
Residual standard error : 35.4 on 68 degrees of freedom
Multiple R-Squared : 0.2077, Adjusted R-squared : 0.1844
F-statistic : 8.914 on 2 and 68 DF, p-value : 0.0003647
Le modèle VAR estimé s’écrit :
Il s’agit d’une version du test de Box et Pierce (1970), corrigée par les degrés de liberté. Ce test
est valable uniquement dans une équation sans autres variables explicatives que des retards de la
variable dépendante. Il permet donc de tester l’autocorrélation du terme d’erreur d’une équation
particulière d’un système VAR. Il est calculé d’après la formule suivante :
s
X rj2
Q∗ (s) = N (N + 1)
j=1
N −j
Le test Q∗ (s) est de type multiplicateur de Lagrange ; il s’agit donc d’un test LM. Sous l’hypothèse
H0 que le terme d’erreur vrai (non observable) de l’équation n’est pas autocorrélé, le test Q∗ (s) est
distribué asymptotiquement selon une Chi-2 à s − p degrés de liberté. On rejette donc l’hypothèse
d’absence d’autocorrélation du terme d’erreur de l’équation si Q∗ (s) prend une valeur réalisée supé-
rieure à la valeur critique d’une Chi-2 à s − p degrés de liberté, au seuil de signification choisi (en
général 5%).
132
b) Test vectoriel portmanteau
Ce test est l’équivalent multivarié du test portmanteau développé pour une seule équation. Le
test vectoriel PB(s) est asymptotique et valable uniquement pour un VAR. Il utilise une correction
pour de petits échantillons. En notant ct le vecteur w ∗ 1 des résidus estimés de toutes les équations
du VAR à la période t, on définit les matrices suivantes, pour tout r, s > 0 :
PN 0
t=1 ct−r ct−r
Crs =
N
où ct−i = 0 pour t − i < 1. Le test vectoriel portmanteau est alors défini de la manière suivante :
s
X tr(C0j C −1 C0j C −1 )
00 00
LB(s) = N 2
j=1
N −j
Il est proposé par Hosking (1980) et Lütkepohl (1991) et est implémenté par exemple dans PcFiml par
Doornik et Hendry ( 1997). Sous l’hypothèse H0 d’absence d’autocorrélation des résidus de toutes les
équations du VAR(p), LB(s) est distribué asymptotiquement selon une Chi-2 à w2 (s − p)/ degrés de
liberté. On rejette donc l’hypothèse d’absence d’autocorrélation des termes d’erreur des w équations
du VAR(p) siLB(s) prend une valeur réalisée supérieure à la valeur critique d’une Chi-2 à w2 (s − p)/
degrés de liberté, au seuil de signification choisi (en général 5%). Le test LB(s) est une version avec
correction du test vectoriel de Box et Pierce qui est :
s
X
P B(s) = N tr(C0j C00−1 C0j C00
−1
)
j=1
# Pour tester l’absence de corrélation sérielle dans les résidus d’un VAR (p), un test de Port-
manteau
ser1 = serial.test(p1, lags.pt = 16, type = "PT.asymptotic")
ser1$serial
Autre tests
Pour tester l’hypothèse d’absence d’autocorrélation du terme d’erreur d’une équation du VAR,
on dispose du test LM habituel d’autocorrélation des résidus. Il repose sur une régression auxiliaire
des résidus sur les variables de droite de l’équation originale et sur les valeurs retardées de ces résidus
à partir d’une période t-r jusqu’à une période t-s. Sous l’hypothèse nulle, le test est distribué selon
une Chi-2 à (s − r + 1) degrés de liberté.
Une valeur trop élevée du test implique donc le rejet de l’hypothèse d’absence d’autocorrélation.
Ce test est valable avec des endogènes retardées, contrairement au test de Durbin et Watson. Il est
possible de généraliser l’approche précédente en testant globalement l’absence d’autocorrélation des
termes d’erreur de toutes les équations du système simultanément.
133
Ce test est basé sur un système auxiliaire, avec des valeurs retardées de résidus à partir d’une
période t r jusqu’à une période t s. Sous l’hypothèse nulle, la version LM de ce test est distribuée
selon une Chi-2 à (s − r + 1)n2 degrés de liberté. Une approximation F du test en question donne
de meilleures propriétés sur de petits échantillons. On peut utiliser un test de White habituel pour
vérifier l’hypothèse d’absence d’hétéroscédasticité du terme d’erreur de chaque équation.
134
est donc permanent et va en s’amortissant.
L’analyse de la réponse impulsionnelle est basée sur la représentation moyenne mobile d’un pro-
cessus VAR (p). Il est utilisé pour étudier les interactions dynamiques entre les variables endogènes.
Les coefficients (i, j) des matrices M sont ainsi interprétés comme la réponse attendue de la variable
yi,t+s à un changement d’unité dans le variable yjt .Ces effets peuvent être cumulés dans le temps
s = 1, 2, . . . . Et donc on obtiendrait l’impact cumulé d’un changement d’unité dans la variable j
à la variable i au temps s.
En dehors de ces coefficients de réponse impulsionnelle, il est souvent concevable d’utiliser des
réponses impulsionnelles orthogonalisées comme alternative. C’est le cas, si les chocs sous-jacents
sont moins susceptibles de se produire isolément, mais plutôt une corrélation contemporaine entre
P
les composantes du processus d’erreur ut , c’est-à-dire les éléments hors-diagonaux de u sont non-
zéro.
yt = εt +1 εt−1 + ...,
Avec εt = P −1 ut pour i = 0, 1, 2, ... .
Incidemment, comme la matrice P est triangulaire inférieure, il s’ensuit que seul un choc dans la
première variable d’un processus VAR (p) exerce une influence sur toutes les autres et que la deuxième
et les suivantes ne peuvent pas avoir un impact direct sur y1t . Il faut tenir compte de ce qui est dit dans
les réponses orthogonales à l’impulsion. Veuillez noter en outre qu’un ordre différent des variables
pourrait produire des résultats différents en ce qui concerne les réponses impulsionnelles.
L’intervalle percentile standard est calculé comme CI s = [S ∗γ , S ∗1−γ ] ou S ∗γ et S ∗1−γ sont les groupes γ2
2 2 2 2
et 1−γ
2
Quantiles des coefficients estimés de réponse impulsionnelle bootstraped Efron & Tibshirani
(1993).
Le choix du sens de l’impact est donc très important et conditionne les valeurs obtenues. Nous
pouvons observer que l’effet d’une innovation s’estompe au cours du temps. Cela caractérise un pro-
cessus VAR stationnaire. Le graphique 1 présente les deux fonctions de réponses impulsionnelles.
Le problème de la corrélation contemporaine des erreurs et donc de l’impact d’un choc sur une va-
riable est traité, d’une manière générale, par la recherche d’une représentation à erreurs orthogonales.
Une fois les chocs initiaux calculés, nous calculons les fonctions de réponse impulsionnelle comme
précédemment, les chocs se répercutent ensuite sur les deux processus en s’amortissant, signe de la
stationnarité du processus VAR.
135
irf.p1 = irf(p1, impulse = "Y1", response = c("Y1","Y2"), boot = T)
irf.p1
plot(irf.p1)
irf.p2 = irf(p1, impulse = "Y1", response = "Y2", boot = T)
plot(irf.p2)
136
- L’amplitude de la cause varie dans le temps.
A partir de ces 3 hypothèses Granger a disposé son test. Il a présenté ce test sur deux variables
X et Y , on dit que Y est causée par X au sens de Granger si l’introduction de valeurs passée de
Xaméliore la qualité statistique de la régression de Y sur ses propres valeurs passés. Autrement dit,Y
est causée par X si les valeurs de Xdonne une information claire sur la valeur actuelle deY tout en
introduisant ses valeurs passées deX, comme des variables explicatives deY . Pour appliquer ce test
on parle de ces deux hypothèses théoriques. Soit le modèle suivant :
La statistique du test suit une loi de Fisher à (c, n − k − 1) degrés de liberté, où : n : Le nombre
de variable explicative.
c : Le nombre de restrictions.
k : Le nombre de variable exogènes
Pour la première équation et pour une erreur α, si on trouve une valeur F supérieur à la valeur
critique au seuil α, donnée par F (c, n − k − 1), on rejette H0 ; c’est-à-dire que l’hypothèse où «X ne
cause pas Y » est rejetée.
Simultanément, si pour la deuxième équation une valeur on trouve une valeur deF inférieur
àF (c, n − k − 1)au seuil α, on ne peut rejeter l’hypothèseH0 , c’est-à-dire que l’hypothèse où « Y ne
cause pas X» ne peut pas être rejetée. Finalement, on obtient que X ne cause pas Y au sens de
Granger et non l’inverse.
137
# Pour tester si Y1(demande) cause de Y2 (prix)
causality(p1, cause = "Y2")
Granger causality H0 : Y2 do not Granger-cause Y1
data : VAR object p1
F-Test = 20.064, df1 = 1, df2 = 136, p-value = 1.572e-05
On rejette H0 cad la demande cause le prix
138
Chapitre 12
La théorie de la cointégration proposée par Granger (1983) et Engel et Granger (1987), est consi-
dérée par beaucoup d’économistes comme un des concepts nouveaux les plus important dans le
domaine de l’économétrie et de l’analyse des séries temporelles. L’analyse de cointégration consiste à
analyser les erreurs entre les composantes non stationnaires de plusieurs séries. Elle permet d’iden-
tifier clairement la relation véritable entre deux variables en recherchant l’existence d’un vecteur de
cointégration et en éliminant son effet.
yt = 5, 92 + 0, 03xt
(8,55) (19,8)
R2 = 0, 94; DW = 0, 057
y = 1 : 30, x = y 2
mod = lm(y ∼ x)
summary(mod)
Les coefficients de régression sont significatifs, la valeur du R2 est élevée, cependant il apparaît évident
que ce modèle a un pouvoir prédictif très mauvais. En effet, sur le plan statistique, la statistique de
Durbin et Watson (proche de 0) présage d’une autocorrélation forte des erreurs1.
Ce premier exemple illustre le danger d’interpréter et d’utiliser une régression entre deux variables
affectées de tendances déterministes de degré différent.
139
12.1.1 Propriétés de l’ordre d’intégration d’une série
Une série est intégrée d’ordre «d» noté xt → I(d) s’il convient de la différence «d» fois afin de la
stationnarité.
- Soit une série x1t stationnaire et une série x2t intègre d’ordre 1
140
Dans le cas général à k variables, on a :
x1t → I(d)
x2t → I(d)
xkt → I(d)
on note Xt = [x1t , x2t , ....., xkt ]
S’il existe un vecteur de cointégration α = [α1t , α2t ...........αkt ] de dimension (k, 1) tel queαxt →
I(d − b). Alors les k variables sont cointégrées et le vecteur de cointégration estα. On note que
xt → CI(d, b) avec b > 0
Notons qu’une grande partie de la littérature sur la cointégration a été consacrée au cas où chaque
variable est intégrée à l’ordre 1. Le principe raison associée est que peu de variables économiques
sont intégrées à un ordre supérieur à 1.
L’idée sous jacente est que si deux variables y et x sont et surcroît cointégration, alors il existe une
combinaison linèaire Z de ces variables qui soit I(0) .Cette combinaison est déterminée par la méthode
de moindres carrée ordinaire.Il s’agit de régression y sur x pour obtenir une estimation de Z. Tester la
cointégration revient à tester la racine unitaire dans Z, pour l’une des méthodes discutées (DF, ADF).
Le test proposé par Engel et Granger présente l’avantage d’être simple à mettre en œuvre, mais
souffre de plusieurs défaillances. En effet, l’estimation de la relation de long terme présuppose qu’une
variable est traitée comme endogène et l’autre comme exogène. Le problème est qu’à priori on peut
régresser aussi bien y sur x.Or selon la normalisation retenue, on peut aboutir à des résultats contra-
dictoires concernant la cointégration des variables. Par exemple, la première régression peut indiquer
que les variables sont cointégrées alors que la seconde peut aboutir à l’absence de la relation de
cointégration. Ceci est d’autant plus gênant que le test de cointégration doit être invariant à la
normalisation retenue. Par ailleurs, cette méthode ne permet pas de tester l’existence de plusieurs
relations de cointégration lorsque le nombre de variable dépasse deux.
141
manière erronée est très élevé. Nous allons tout d’abord examiner le cas à deux variables : test de
cointégration et estimation du modèle à correction d’erreur.
yt = α̂ + β̂xt + et
ECM
Étape 2 : estimation par les MCO de la relation du modèle dynamique (court terme) :
142
∆yt = α1 ∆xt + α2 et−1 + µt ; α2 < 0
Démonstration
yt = α0 + α1 yt−1 + α2 xt + α3 xt−1 + vt
yt = axt + b + εt
Car à long terme, on a yt−1 = yt ,xt−1 = xt et la dynamique de court terme devient à long terme :
yt = α0 + α1 yt + α2 xt + α3 xt + vt
D’où (1 − α1 )y t = (α2 + α3 )xt + α0 + vt
yt = axt + b + εt
Avec = α1−α
2 +α3
1
α0
, b = 1−α 1
vt
,εt = 1−α 1
.
Le MCE s’obtient à partir de la dynamique de court terme :
yt = α0 + α1 yt−1 + α2 xt + α3 xt−1 + vt
yt − yt−1 = α0 + α1 yt−1 − yt−1 + α2 xt − α2 xt−1 + α2 xt−1 + α3 xt−1 + vt
4y t = (α1 − 1) yt−1 + α2 (xt − xt−1 ) + α0 + (α2 + α3 )xt−1 + vt
4y t = − (1 − α1 ) yt−1 + α2 (xt − xt−1 ) + α0 + (α2 + α3 )xt−1 + vt
α2 + α3 α0
4y t = − (1 − α1 ) yt−1 − xt−1 − + α2 4xt + vt
1 − α1 1 − α1
143
12.4.1 Exemple
Soit deux séries statistiques yt et xt dont les observations sont présentées au tableau 1. On
demande d’estimer la relation entre ces deux variables (yt = a0 + a1 xt + et ) en testant une éventuelle
cointégration (dans ce cas estimer le modèle à correction d’erreur).
La première étape consiste à étudier les propriétés des deux séries en termes de stationnarité.
La comparaison des t calculés aux t lus (Tables de MacKinnon) aux valeurs critiques indique que les
deux séries yt et xt sont non stationnaires en niveau1. Des tests similaires sur les différentes premières
de yt et xt indiquent qu’elles sont stationnaires, les deux séries étudiées sont I (1), il existe donc un
risque de cointégration.
Le test de cointégration est effectué à partir du résidu d’estimation du modèle :
Soit : yt = 0.55xt + 10, 38 + et
(6,3) (41,46)
2
n = 30; R = 0, 58; (.) = tdeStudent
Nous pouvons vérifier que le résidu est bien stationnaire, il existe donc un risque de cointégration
entre les deux variables.
Nous procédons donc maintenant à l’estimation du modèle à correction d’erreur.
Nous calculons, d’abord, le résidu (provenant du modèle précédent) décalé d’une période, soit :et−1 =
yt−1 − 0, 55xt−1 − 10, 38 Puis nous estimons (par les MCO) le modèle :∆yt = α1 ∆xt + α2 et−1 + µt
n = 29; R2 = 0, 60; (.) = tdeStudent. Le coefficient (terme de rappel) de et−1 est bien significative-
ment négatif, la représentation à correction d’erreur est validée.
144
stationnaire elles sont alors cointégrées. L’estimation par lesM COdu modèle permet de calculer le
résidu.
et = e1t + e2t = yt − α̂0 − α̂1 x1t + x2t − γ̂0 − γ̂1 x3t → I(0)
Nous obtenons, dans ce cas de figure, un autre vecteur de cointégration possible :
[1, −α̂0 − γ̂0 , −α̂1 1, −γ̂1 ]
D’une manière générale, dans un modèle à une variance à expliquer et k variables explicatives,
il peut exister k vecteurs de cointégration linéairement indépendants. Le nombre de vecteurs de
cointégration linéairement indépendants est appelé le rang de la cointégration.Si les variables sont
de même ordre d’intégration – I (1) par exemple – l’existence d’un seul vecteur de cointégration
est possible ; en revanche, si les séries ne sont pas toutes intégrées du même ordre, nous pouvons
être certains que le vecteur de cointégration n’est pas unique. De manière pratique, pour tester une
éventuelle cointégration entre plusieurs variables, il convient tout d’abord de la tester sur l’ensemble
des k + 1 variables, puis – en cas de cointégration – de la tester par combinatoire entre les variables
145
Le coefficient γt doit être, comme dans le cas à une variable explicative, significativement négatif.
Ce pendant, le plus souvent, le vecteur de cointégration n’est pas unique et la méthode d’Engel-
Granger n’est plus validée. En effet, les estimateurs des MCO ne sont plus consistants quelque soient
les vecteurs de cointégration. Nous devons, dans ce cas, faire appel à la représentation vectorielle à
correction d’erreur (VECM, « Vector Error Correction Model »).
Yt = A0 + A1 Yt−1 + A2 Yt−2 +
avec :
Yt : vecteur de dimension (k ,1) constitué des k variables (y1t , y2t ..., ykt ), A0 : vecteur de dimension
(k , 1),
Ai : matrice de dimension (k , k).
Ce modèle peut s’écrire en différence première :
Yt − Yt−1 = A0 + A1 Yt−1 − A1 Yt−1 + A2 Yt−2 +
∆Yt = A0 + (A1 − I)Yt−1 + A2 Yt−2 +
Afin de faire apparaître des différences premières à droite de l’équation, nous ajoutons et soustrayons
A1 Yt−2 − Yt−2 de la manière suivante
∆Yt = A0 + (A1 − I)Yt−1 + A1 Yt−2 − Yt−2 − A1 Yt−2 + Yt−2 + A2 Yt−2 +
∆Yt = A0 + A1 Yt−1 − Yt−1 + A1 Yt−2 − Yt−2 − A1 Yt−2 + Yt−2 + A2 Yt−2 +
avec :
Yt : vecteur de dimension (k ,1) constitué des k variables (y1t , y2t ..., ykt ), A0 : vecteur de dimension
(k , 1),
Ai : matrice de dimension (k , k).
146
Ce modèle peut s’écrire en différences premières de deux manières :
∆Yt = A0 + (A1 − I)∆Yt−1 + (A1 + A2 − I)∆Yt−2 + ... + (Ap−1 + .... + A2 + A1 − I)∆Yt−p+1 + πYt−p +
ou encore en fonction de Yt−1 :
∆y1,t = a10 + b11 ∆y1,t−1 + b12 ∆y2,t−1 + b13 ∆Y3,t−1 + α1 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 ) + 1t
∆y2,t = a20 + b21 ∆y1,t−1 + b22 ∆y2,t−1 + b23 ∆Y3,t−1 + α2 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 ) + 2t
∆y3,t = a30 + b31 ∆y1,t−1 + b32 ∆y2,t−1 + b33 ∆Y3,t−1 + α3 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 ) + 3t
Nous avons normalisé l’équation de cointégration par rapport au coefficientde y1,t .
Selon les caractéristiques des données traitées, ce processus VAR peut connaître les variantes sui-
vantes :
– existence d’une constante soit dans la relation de cointégration (β0 ) soit dans le VAR (ai0 ), – exis-
tence d’une tendance (t = 1,2,. . . ,n) dans le VAR et/ou dans la relation de cointégration,
– ou encore intégrer des variables exogènes de type indicatrice pour corriger d’un mouvement saison-
nier.
Sur ce type de modèle, nous ne pouvons pas appliquer la méthode des MCO car nous avons des pro-
blèmes d’identification similaires à ceux évoqués lors du chapitre 8 concernant les modèles à équations
simultanées. Il convient d’employer une méthode du maximum de vraisemblance1.
147
Étape 1 : calcul de deux résidus ut et vt
Nous effectuons deux régressions :
Première régression : ∆Yt = A0 + A1 ∆Yt−1 + A2 ∆Yt−2 + ... + Ap ∆Yt−p + ut
Deuxième régression : Yt = A0 + A1 ∆Yt−1 + A2 ∆Yt−2 + ... + Ap ∆Yt−p + vt
Avec Yt = [y1,t y2,t , ....., yk,t ]
Nous avons les mêmes variables explicatives, seule la spécification du bloc de la variable à expliquer
est modifiée. ut et vt sont donc les matrices des résidus de dimension (k, n) avec k = nombre de
variables, n = nombre d’observations.
Étape 2 : calcul de la matrice permettant le calcul des valeurs propres
Nous calculons quatre matrices des variances-covariances de dimension (k, k) à partir des résidus ut
et vt .
n n n n
1X 0 1X 0 1X 0 1X 0
Σuu = ut ut ; Σvv = vt vt ; Σuv = ut vt ; Σvu = vt ut
n t=1 n t=1 n t=1 n t=1
Puis nous extrayons les k valeurs propres de la matrice M de dimension k, k calculée de la manière
suivante :
M = Σ−1 −1 −1
vv Σuv Σuu
k
X
λtrace = −n Ln(1 − λi )
i=r+1
148
Pour mener ce test, Johansen propose cinq spécifications concernant soit les vecteurs cointégrants
soit les séries (le VAR proprement dit) :
Absence de tendance linéaire dans les données (les séries sont toutes DS sans dérive) : a) Absence
d’une tendance linéaire dans les séries et d’une constante dans les relations de cointégration (la
constante dans la relation de long terme est non significative).
∆y1,t = b11 ∆y1,t−1 + b12 ∆y2,t−1 + b13 ∆Y3,t−1 + α1 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 ) + 1t
b) Absence d’une tendance linéaire dans les séries mais présence d’une constante dans les relations
de cointégration (la constante dans la relation de long terme est significative).
∆y1,t = b11 ∆y1,t−1 + b12 ∆y2,t−1 + b13 ∆Y3,t−1 + α1 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 ) + 1t
Présence d’une tendance linéaire dans les données (au moins une série est un DS avec dérive) :
c) Présence d’une tendance linéaire dans les séries et d’une constante dans les relations de cointégra-
tion.
∆y1,t = a10 + b11 ∆y1,t−1 + b12 ∆y2,t−1 + b13 ∆Y3,t−1 + α1 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 ) + 1t
d) Présence d’une tendance linéaire dans les séries et dans les relations de cointégration (au moins
une série est un TS).
∆y1,t = a10 + b11 ∆y1,t−1 + b12 ∆y2,t−1 + b13 ∆Y3,t−1 + α1 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 + ct) + 1t
Présence d’une tendance quadratique dans les données :
e) Présence d’une tendance quadratique dans les séries et d’une tendance linéaire dans les relations
de cointégration
∆y1,t = a10 + bt + b11 ∆y1,t−1 + b12 ∆y2,t−1 + b13 ∆Y3,t−1 + α1 (y1,t−1 − β2 y2,t−1 − β3 y3,t−1 − β0 + ct) + 1t
Le choix d’une de ces spécifications s’effectue en fonction des données et de la forme supposée de la
tendance (une analyse des propriétés stochastiques des séries ou un simple examen visuel des gra-
phiques permettent de se déterminer).
λmax = −nLog(1 − λ + 1) r = 0, 1, 2. . .
Le test s’effectue comme précédemment de manière séquentielle par exclusion d’hypothèses alterna-
tives. En cas de divergence des deux tests (valeur propre maximum et trace), nous privilégions le
test de la trace dont la puissance est la plus élevée.
149
12.7 Synthèse de la procédure d’estimation
Nous essayons ici de synthétiser les grandes étapes relatives à l’estimation d’un modèle VECM.
Étape 1 : Détermination du nombre de retards p du modèle (en niveau ou en Log) selon les critères
AIC ou SC .
Étape 2 : Estimation de la matrice π et test de Johansen permettant de connaître le nombre de
relations de cointégration (les logiciels proposent un certain nombre de spécifications alternatives,
telles que l’existence d’un terme constant dans la relation de cointégration, contraindre A0 = 0 ,
l’existence d’une tendance déterministe, etc.).
Étape 3 : Identification des relations de cointégration, c’est-à-dire des relations de long terme entre
les variables.
Étape 4 : Estimation par la méthode du maximum de vraisemblance du modèle vectoriel à correction
d’erreur et validation à l’aide des tests usuels : significativité des coefficients et vérification que les
résidus sont des bruits blancs (test de Ljung-Box), tests d’exogénéité faible.
Enfin, nous pouvons vérifier que l’estimation par les MCO de la relation de long terme fournit des
résultats à peu près similaires (en termes de significativité et de valeurs estimées des coefficients) à
ceux obtenus par la méthode du maximum de vraisemblance.
150
Test type : trace statistic , with linear trend
Eigenvalues (lambda) :
6.048611e-01 2.569525e-01 1.390632e-01 -4.180306e-17
Les trois valeurs propres de la matrice π, estimée par le maximum de vraisemblance, sont égales à
λ1 = 0, 605; λ2 = 0, 257; λ3 = 0.139.
Premier test : Rang de la matrice π égal 0 (r = 0), soit H0 : r = 0 contre H1 : r > 0.
Values of test statistic and critical values of test :
test 10pct 5pct 1pct
r <= 2 | 4.19 7.52 9.24 12.97
r <= 1 | 12.51 17.85 19.96 24.60
r = 0 | 38.51 32.00 34.91 41.07
Calculons la statistique de Johansen :λtrace = −n ki=r+1 Ln(1 − λi )pourr = 0
P
151
fication est rejetée du fait que les trois constantes des trois équations ne sont pas significativement
différentes de 0. L’estimation finale du VECM sur 28 observations est donc la suivante :
vecm.r1 = cajorls(sjf.vecm1, r = 1)
vecm.r1
lm(formula = substitute(form1), data = data.mat)
Coefficients :
Y1.d Y2.d Y3.d
ect1 -0.85 0.32 0.0089
Y1.dl1 -0.08 0.43 0.54
Y2.dl1 0.02 -0.62 -0.46
Y3.dl1 0.20 -0.25 -0.25
$beta
ect1
Y1.l2 1.00
Y2.l2 -0.87
Y3.l2 -0.54
const 12.92
152
Bibliographie
[1] Bourbonnais, R ; Terraza,M (2008), Analyse des séries temporelles, Dunod, 2e édition
[2] Ihaka, R ; Gentleman, R. (1996) R : A language for Data Analysis and Graphics. Journal of
Computational and Graphical Statistics , 5(3), 299–314.
[3] Cribari-Neto, F. Zarkos, S. G. (1999), R : yet another econometric programming environment, ,
Journal of Applied Econometrics 14(3), 319
[4] Cryer et Chan (2008) Time Series Analysis with Applications in R, Springer
[5] Mignon.V et S. Lardic (2002), Analyse des séries temporelles, Econometrica
[6] Pirotte. A et G. Bresson (1995) Économétrie des séries temporelles, PUF
[7] Pfaff B (2006) Analysis of Integrated and Cointegrated Time Series with R, Springer
[8] Fox, J.(2017) Using the R Commander : A Point-and-Click Interface for R, 1st ed. ; Chapman
and Hall/CRC Press : BocaRaton, FL.
[9] Mestiri, S. (2019) (2019) How to use the R software, MPRA Paper, University Library of
Munich, Germany.
[10] Mestiri, S. (2023) Initiation à l’utilisation du logiciel R Éditions universitaires européennes
[11] Mestiri,S. (2021) Simulation des données en utilisant logiciel R. Available at HAL Working
Papers hal-03448651
[12] Mestiri, S. (2024) Machine Learning Techniques in Financial Applications. Journal of
Research, Innovation and Technologies, Vol 3, 1(5), 5-24.
[13] Mestiri, S. (2024) ARDL modeling using R software. Journal of Current Trends in
Computer Science Research, Vol. 3 No 1, pages, 01-05
[14] Venables, W. Ripley, B. (1999), Modern applied statistics with S-Plus, 3rd edn, SpringerVerlag,
New York.
[15] Shumway R et Stoffer D. (2006) Time Series analysis and its applications with R examples,
springer
[16] Greene.W (2005), Econometrics, Pearson Education
153