Rappel sur les séries chronologiques univariées

Série temporelle →suite d’observations yt , t  T  d’une variable y à différentes dates t

L’indice temps→ l’année, le trimestre l’heure, le mois, le jour, …
Exemples → le PIB d’un pays, les ventes d’une entreprise, le cours d’une action…
Processus stochastique →suite ordonnée de variables aléatoires yt (),   , t  T 
Une série temporelle est une réalisation d’un processus aléatoire
Une étude d’une série temporelle consiste à
 analyser, de décrire et d’expliquer un phénomène
( Quel processus stochastique qui a généré la série observée)
 prévoir ce phénomène dans le futur à partir des observations
 tirer des conséquences pour des prises de décision
Représentation graphique d’une série temporelle→Pour une étude descriptive, la
représentation graphique d’une série temporelle est nécessaire : yt en ordonnées et t en abscisses.
Exemple : le nombre de passagers (en milliers) par mois dans les transports aériens, de 1949 à 1960
(disponible sous R :AirPassengers)
ts.plot(AirPassengers,xlab="t",ylab="y")
600
500
400
y

300
200
100

1950 1952 1954 1956 1958 1960

Indices descriptifs d’une série temporelle qui ‘‘résument’’ la série yt tn1

1 n
Indices de tendances centrales→ La moyenne y   yt et la médiane.
n t 1

1 n 2
Indices de dispersion→ La variance empirique  ( yt  y ) ,
n t 1

Indices de dépendance → l’auto-covariance et l’auto-corrélation empiriques →

n
 ( yt  y )( yt  h  y )
Ch  t  h 1
C
et rh  h (C0 est la variance empirique)
nh C0

1
 n n
(1 / n)  ( yt  y )3 (1 / n)  ( yt  y ) 4
Indices d’asymétrie et d’aplatissement → ˆ1  t 1 et ˆ2  t 1 3
C03 / 2 C02

Composantes principales d’une série temporelle

Tendance → ft représente les effets permanents
(Exemple : yt  1  0.05t   t ,  t ~ N (0,1) et n  250)

15
10
y

5
0

0 50 100 150 200 250

Saisonnalité → st représente les effets périodiques

 *t
(Exemple : yt  2 * sin(  4)   t ,  t ~ N (0,1) et n  100
6
3
2
1
0
y

-1
-2
-3
-4

0 20 40 60 80 100

Une composante résiduelle →  t représente les effets imprévisibles

Un modèle de série temporelle→ est une équation reliant ces composantes. Parmi ces modèles :
- le modèle additif : yt  f t  st   t (Les droites reliant les minimax et les maximas sont parallèles)
-le modèle multiplicatif : yt  f t  st   t (cas contraire)
Détection, estimation et élimination de la tendance
Détection : - le graphe ; un test
Estimation :
- Estimation paramétrique : yt  f t   t , ft  a  bt
- Estimation non paramétrique : Méthode des moyennes mobiles

2
Elimination:
- soustraire la tendance à la série initiale zt  yt  ft
- recourir aux différences d’ordre 1, zt  yt  yt 1
Détection, estimation et élimination de la saisonnalité
Détection : le graphe ; calcul du coefficient d’autocorrélation empirique rh ; un test
Estimation :
- Estimation paramétrique : yt  st   t , st =sin (ωt + ϕ)
- Estimation non paramétrique : → calcul des coefficients saisonniers corrigés
Elimination :
 pour un modèle additif, ydsi, j  yi, j  s j ( ydsi, j  yi, j / s j  cas multiplicatif)

 ydst  yt  yt  p ( application de filtres différences de type p)

Appliquer deux filtres différences si l’on veut, à la fois, éliminer la tendance et la saisonnalité
Prévision→ Soit une série y1 ,..., yn . Prévoir yn  h , non encore observée.
- la prévision est notée yˆ n (h) ,
- l’erreur de prévision est donnée par la différence en  h  yˆ n (h)  yn  h
Méthodes de prévision
- Modèles déterministes : Exemple : yt  a  bt   t , la prévision de yt 1 sera aˆ  bˆ(t  1)
- Les lissages exponentiels
- Les modèles de type ARMA
Concepts fondamentaux d’une série temporelle
Processus stochastique et série temporelle
(, A, P) un espace de probabilité et (Yt ) un processus stochastique (T  Z)

- pour tout t  T , Yt est une variable aléatoire de  dans R

- pour tout    , (Yt ( )  yt , t  T ) est appelée série temporelle

 E ( yt )    Cste
 2
Stationnarité→ ( yt ) est faiblement stationnaire →  V ( yt )     0  
cov( y , y
 t t  h )   h , t , h

( yt ) est fortement stationnaire→pour tout t1  t 2  ...  t n (ti  Z ) et h  Z

les suites ( yt1 h ,..., ytn  h ) et ( yt1 ,..., yt n ) ont même distribution

Bruit blanc : Soit une suite ( t ) de variables aléatoires

- Variables  t centrées, identiquement distribuées et indépendantes → ( t ) est un bruit blanc fort

- Variables  t centrées, identiquement distribuées et non corrélées → ( t ) est un bruit blanc faible

E( t )  0,V ( t )   2 et  (h)  cov( t ,  t  h )  0 pour h  0 (Si  t ~ N (0,  ) alors ( t ) est un BB gaussien)

3
  ( h)
Fonction d’autocorrélation : h  → corrélation de la série avec elle-même décalée de h
 (0)

périodes.
Fonction d’autocorrélation partielle → corrélation de yt et yt h
l’influence des variables yt 1 , yt 2 …, yt h1 ayant été retirée.
Test de corrélation : H 0 :  h =0 contre H1 :  h  0
n2
Pas de corrélation si t  tn2 ( / 2) où, t  rh
1 r 2
h
Test de bruit blanc ( Box-Pierce) : H 0 : 1  ...   h  0
h
Nous acceptons l’hypothèse de bruit blanc H 0 si Q   (2h) ( ) , où Q  n  rk2 , (ou p-value > α)
k 1
Test d’homoscédasticité (Test ARCH) . Soit et2  a0  a1et21  ...  aq et2 q  t ...(1)

L’hypothèse d’homoscédasticité H 0 : a1  ...  aq  0 est acceptée si nR 2   (2q) ( ) , (ou p-value > α)

où n est le nombre d’observations et R2 est le coefficient de détermination du modèle (1).

Test de normalité( Jarque et Bera) : H 0 : 1  0 et  2  0

Nous acceptons l’hypothèse de normalité si JQ   (22) ( ) (ou p-value > α)

ˆ  0 ˆ  0
où , JQ  g12  g 22 avec g1  1 et g 2  2
6/ n 24 / n

Les processus AR, MA , ARMA

Processus autorégressif AR(p) → yt  c  1yt 1  2 yt 2  ...  p yt  p   t ,où  t ~ BB(0, 2 )

Processus moyenne mobile d’ordre MA(q)→ yt  m   t  1 t 1  2 t 2  ... q t q ,

Processus autorégressif moyenne mobile ARMA(p,q)

yt  c  1yt 1  2 yt 2  ...  p yt  p   t  1 t 1   2 t 2  ...  q t q

AR(p) → f.a.p s’annule à partir du rang p  1

MA(q) →f.a s’annule à partir du rang q+1
Identification, estimation, validation et prévision
Box et Jenkins → ont proposé une procédure en plusieurs étapes
- Transformer yt afin d’éliminer d’éventuelles non-stationnarité
Par exemple, si la variance de la série croit avec t on peut utiliser une transformation de type
logarithmique ou si la série présente une saisonnalité on retire la composante saisonnière avant tout
traitement statistique et l’ajouter à la fin du traitement
- Identifier le modèle ARMA(p,q) sur la base des caractéristiques de yt ( autocorrélations simple et
partielle)

4
- Estimer les paramètres du modèle c , i (i  1,..., p) ,  j ( j  1,...,q) et  2 d’un modèle ARMA à partir

d’observations y1,..., yn , ( Méthodes : exemple → MV)

-Tester la validité et l’adéquation du modèle aux données,
Des tests sont appliqués sur les paramètres (student) et sur les résidus ( BBG)
Absence d’autocorrélation → Test de Box Pierce
Test d’homoscédasticité→Test ARCH :
Normalité →tests de Jarque-Bera, Shapiro-Wilk, …
- Re-spécifier le modèle si besoin ou faire de la prévision.
Cas où plusieurs modèles valides→ meilleur modèle selon critère AIC ou performance en prévision
La prévision, yˆt  h = E( yt h It ) , où I t est l’ensemble d’information disponible à la date t

La non stationnarité. Tests de stationnarité.

La plupart des séries économiques et financières → non stationnaires

La non stationnarité : déterministe ou stochastique
yt  f t  zt , → une non stationnarité déterministe TS où ft dépend de t et zt stationnaire peut être

modélisée par un ARMA→ stationnariser par régression sur le temps.

zt  yt  yt 1 est stationnaire→une non stationnarité stochastique DS ( yt ~ I (1) :Processus intégré

d’ordre un) → stationnariser par passage aux différences

Généralement, ( yt ) est un processus intégré d’ordre d ( yt ~ I (d )) si (1  L) d yt est stationnaire

Les exemples les plus connus : marche aléatoire sans dérive ( yt  yt 1  zt ) et marche aléatoire avec
dérive ( yt  c  yt 1  zt ) → Présence de racine (s) unitaire (s) pour non stationnairité stochastique
Il existe plusieurs tests de racine unitaire :Test de Dickey-Fuller simple, Test de Dickey-Fuller
augmenté, Test de Phillips et Perron, Test KPSS,….
- Le test de Dickey-Fuller simple :Trois modèles de base sont considérés,
( M 1 ) : yt  yt 1   t …(1), ( M 2 ) : yt  yt 1     t , …(2), ( M 3 ) : yt  yt 1    t   t …(3)
Dans le modèle simple (1), le test H0 :   1 ( présence de racine unitaire) contre H1 :   1 (modèle
ˆ ˆ
AR(1) stationnaire) → accepter H 0 si tˆ  ttabulé où t    1 (où t  , si     1)
ˆ (ˆ) ˆ ( ˆ )
La stratégie du test de Dicker-Fuller est donné par un tableau résumant toutes les règles de décision à
prendre. A cause du rejet du bruit blanc, une extension du test précédent permet d’intégrer la présence
d’autocorrélétion. On applique le test de Dickey-Fuller augmenté dans les trois modèles suivants :
p p
( M 1 ) : yt  yt 1    j yt  j  t …(1), ( M 2 ) : yt    yt 1    j yt  j  t , … (2),
j 1 j 1
p
( M 3 ) : yt    t  yt 1    j yt  j  t …(3)
j 1

5
Processus ARIMA et SARIMA

- ( yt ) est un processus ARIMA(p,d,q) si zt  (1  L) d yt est un processus ARMA(p,q)

- ( yt ) sera un processus SARIMA(p,d,q)  (P,D,Q) avec période s si zt  (1  L) d (1  Ls ) D yt est un
processus ARMA,où d différence, s saison, D différence saisonnière, l’opérateur de différentiation
(1  Ls ) D élimine la saisonnalité l’opérateur (1  L) d élimine les tendances.

Exemple empirique d’un modèle ARMA

Considérons la série Wind, correspondante à des données journalières du vitesse du vent du

mois de mai jusqu’au mois de septembre, obtenue à partir de data (airquality) où Wind= airquality[,3].
La série d’étude, Wind1 est la série Wind privée de ses 40 premières valeurs : Wind1=Wind[41 :153].
1. Représenter graphiquement la série Wind1 et donner des commentaires concernant ses
composantes.
2. Tracer graphiquement les fonctions d’auto-corrélation simple et partielle et donner des
commentaires concernant la stationnarité de la série et la spécification du modèle.
3. Ajuster un modèle ARMA pour la série Wind1.
4. Utiliser les tests de Box-Pierce, d’Arch et de Jarque-Bera pour la série des résidus obtenue.
5. Calculer les prévisions pour un horizon h=5.
****************************************************************
data(airquality) ; Wind=airquality[,3] ; Wind1=Wind[41 :153]
1. Stationnarité de la série
Graphe de la série : ts.plot(Wind1)
20
15
Wind1

10
5

0 20 40 60 80 100

Time

Fig1
La série parait stationnaire. Confirmons la stationnarité de la série Wind1par le test de Dickey-Fuller.
Commencer par télécharger le package urca .
library(urca);
Considérons le modèle (1) avec trend ( tendance).
x1=ur.df(Wind1,type='trend'); summary(x1)
Augmented Dickey-Fuller Test Unit Root Test # Test regression trend
Coefficients: Estimate Std. Error t value Pr(>|t|)

6
 tt 0.004857 0.009623 0.505 0.615
Value of test-statistic is: -5.1539 8.8675 13.301
Critical values for test statistics: 1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
tc =-5.15 < tau3= -3.99→accepter H1: rau<0
p-value de beta= 0.61> 0.05→accepter H0 :beta=0
On passe au modèle (2) avec drift (constante)
x2=ur.df(Wind1,type='drift'); summary(x2)
Augmented Dickey-Fuller Test Unit Root Test # Test regression drift
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.28548 1.07253 4.928 3.02e-06 ***
Value of test-statistic is: -5.1508 13.2655
Critical values for test statistics: 1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
tc =-5.15<tau2= -3.46→accepter H1: rau<0
p-value de alpha= 3.02*10-6 < 0.05→rejeter H0 :alpha=0
On s’arrête au modèle (2). La série Wind1 est donc stationnaire et est de la forme I(0)+c.
2. Identification
Fonction d’autocorrélation simple et partielle :
acf(Wind1,25) ; pacf(Wind1,25) ;
Series Wind1 Series Wind1
1.0

0.3
0.8

0.2
0.6

Partial ACF

0.1
ACF

0.4

0.0
0.2

-0.1
0.0
-0.2

-0.2

0 5 10 15 20 25 5 10 15 20 25

Lag Lag

Fig2 Fig3
Nous prenons pmax=1 et qmax=2 et donc cinq modèles possibles pour la série Wind1 à
savoir :AR(1), MA(1), ARM(1,1), MA(2) et ARMA(1,2).
3. Estimation
Fit1=armaFit(~arma(1,0),data=Wind1); summary(Fit1)
Coefficient(s): Estimate Std. Error t value Pr(>|t|)
ar1 0.33999 0.08818 3.856 0.000115 ***
intercept 9.45493 0.45417 20.818 < 2e-16 ***
sigma^2 estimated as: 10.24; log likelihood: -291.86; AIC Criterion: 589.71

7
Fit2=armaFit(~arma(0,1),data=Wind1); summary(Fit2)
Coefficient(s): Estimate Std. Error t value Pr(>|t|)
ma1 0.26377 0.07999 3.298 0.000975 ***
intercept 9.44643 0.38603 24.471 < 2e-16 ***
sigma^2 estimated as: 10.58; log likelihood: -293.67;AIC Criterion: 593.33
Fit3=armaFit(~arma(1,1),data=Wind1); summary(Fit3)
Coefficient(s): Estimate Std. Error t value Pr(>|t|)
ar1 0.5859 0.1563 3.749 0.000178 ***
ma1 -0.2735 0.1754 -1.559 0.118990
intercept 9.4661 0.5194 18.226 < 2e-16 ***
sigma^2 estimated as: 10.08; log likelihood: -290.94; AIC Criterion: 589.88
Fit4=armaFit(~arma(0,2),data=Wind1); summary(Fit4)
Coefficient(s): Estimate Std. Error t value Pr(>|t|)
ma1 0.24259 0.09814 2.472 0.0134 *
ma2 0.19654 0.09893 1.987 0.0470 *
intercept 9.44784 0.43076 21.933 <2e-16 ***
sigma^2 estimated as: 10.2; log likelihood: -291.63; AIC Criterion: 591.26
Fit5=armaFit(~arma(1,2),data=Wind1); summary(Fit5)
Coefficient(s): Estimate Std. Error t value Pr(>|t|)
ar1 0.4469 0.2048 2.182 0.0291 *
ma1 -0.1788 0.2015 -0.887 0.3749
ma2 0.1635 0.1231 1.329 0.1839
intercept 9.4610 0.5237 18.067 <2e-16 ***
sigma^2 estimated as: 9.938; log likelihood: -290.19;AIC Criterion: 590.37
4. Validation
Première étape : Test de student sur les paramètres
Selon le test de signification sur les paramètres, les modèles 3 et 5 sont rejetés.
Deuxième étape : Tests sur les résidus
Nous testons les résidus obtenus des modèles non rejetés pour savoir s’ils constituent un bruit blanc
gaussien. d’autocorréllation ).
res1=residuals(Fit1) ;ts.plot(res1) ; acf(res1)
Series res1
10

1.0
0.8
5

0.6
ACF
res1

0.4
0

0.2
0.0
-5

-0.2

0 20 40 60 80 100 0 5 10 15 20 25

Time Lag

Fig4 Fig5

8
Test d’absence d’autocorréllation ( Nous utilisons le test de Box et Pierce)
Box.test(res1,25)
X-squared = 21.9676, df = 25, p-value = 0.6376
p-value > 0.01 →hypothèse d’autocorréllation rejetée
Test d’homoscédasticité des résidus (Nous utilisons le test ARCH)
library(FinTS); ArchTest(res1);
Chi-squared = 21.104, df = 12, p-value = 0.04888
p-value > 0.01→ hypothèse d’homoscédasticité acceptée
Test de normalité des résidus ( Nous utilisons le test de Jarque et Bera)
library(tseries) ;
jarque.bera.test(res1)
X-squared = 0.9232, df = 2, p-value = 0.6303
p-value > 0.01→ hypothèse de normalité acceptée
Les résidus vérifient bien un bruit blanc gaussien pour α=0.01.
D’une manière similaire, on peut vérifier que les résidus obtenus des estimations selon les modèles
MA(1) et MA(2) suivent aussi un bruit blanc gaussien.
Comparaison des trois modèles candidats selon le critère AIC
p 1 0 0
q 0 1 2
AIC 589.71 593.33 591.26

Nous choisissons le modèle AR(1) car il possède une valeur AIC minimum.
On peut utiliser l’instruction auto.arima du package forecast pour choisir le modèle optimal.
library(forecast) ;auto.arima(Wind1)
Series: Wind1 ; ARIMA(1,0,0) with non-zero mean
Coefficients: ar1 intercept
0.3400 9.4549
s.e. 0.0882 0.4542
sigma^2 estimated as 10.24: log likelihood=-291.86; AIC=589.71

5. prevision
predict(Fit1,5)
Low 95 Low 80 Forecast High 80 High 95
114 3.8773 6.0486 10.1502 14.2519 16.4231
115 3.0658 5.3591 9.6913 14.0235 16.3169
116 2.8702 5.1772 9.5353 13.8934 16.2004
117 2.8126 5.1212 9.4823 13.8433 16.1519
118 2.7940 5.1028 9.4642 13.8256 16.1344

9
 ARIMA(1,0,0) with method: CSS-ML

16
14
12
Series: x

10
8
6
4
2

70 80 90 100 110 120

Time

Fig6

Exemple empirique d’un modèle ARIMA

Nous nous intéressons à la modélisation du logarithme du PNB américain sur la période 1970-
1987 en dollars. Les observations sont obtenues à partir de data(USeconomic) du package tseries
(Utiliser les instructions suivantes, LPNB = USeconomic[,2]; LPNB1=LPNB[65:136]).
1. Représenter graphiquement la série LPNB1et sa fonction d’autocorrélation pour 20 retards.
Donner des commentaires.
2. Tester la stationnarité de la série différenciée d’ordre 1 dLPNB1 par le test Dickey-Fuller
(dLPNB1=diff(LPNB1).
3. Ajuster un modèle ARIMA pour la série LPNB1.
4. Utiliser les tests de Box-Pierce, Arch et Jarque-Bera pour la série des résidus
obtenus. Que peut-on déduire ?
5. Effectuer une prévision pour l’année 1988.
****************************************************************
library(tseries) ;
data(USeconomic); LPNB = USeconomic[,2]; LPNB1=LPNB[65:136];
1. Le graphe de la série LPNB1 et sa fac pour 20 retards
ts.plot(LPNB1) ;acf(LPNB1,20)

10
 Series LPNB1

1.0
8.2

0.8
0.6
8.1
LPNB1

ACF

0.4
8.0

0.2
7.9

0.0
-0.2
7.8

0 10 20 30 40 50 60 70 0 5 10 15

Time Lag

La présence d’une tendance et la non décroissance rapide de la fonction d’autocorrélation indiquent la

non stationnarité de la série LPNB1. Nous pouvons confirmer ce résultat par le test de Dickey-Fuller )
2. Test de stationnarité de la série différenciée dLPNB1
dLPNB1=diff(LPNB1) ;pacf(dLPNB1,20)

Series dLPNB1
0.2
0.1
Partial ACF

0.0
-0.1
-0.2

5 10 15

Lag

Utilisons le test de Dickey-Fuller pour la série différenciée

x1=ur.df(dLPNB1,lag=1,type='trend'); summary(x1)
# Augmented Dickey-Fuller Test Unit Root Test # Test regression trend
Coefficients: Estimate Std. Error t value Pr(>|t|)
tt 1.778e-05 6.492e-05 0.274 0.785094
Value of test-statistic is: -4.0967 5.6138 8.4163
Critical values for test statistics: 1pct 5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2 6.50 4.88 4.16
phi3 8.73 6.49 5.47
tc =-4.09 <tau3= -4.04→accepter H1: rau<0
p-value de beta= 0.78> 0.05→accepter H0 :beta=0
On passe au modèle (2)

11
 x2=ur.df(dLPNB1,lag=1,type='drift'); summary(x2)
# Augmented Dickey-Fuller Test Unit Root Test # Test regression drift
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.004268 0.001631 2.617 0.010986 *
Value of test-statistic is: -4.1226 8.5024
Critical values for test statistics: 1pct 5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1 6.70 4.71 3.86
tc =-4.12<tau2=-3.51→accepter H1: rau<0
p-value de alpha= 0.01 < 0.05→rejeter H0 :alpha=0→ dLPNB1 est de la forme I(0)+c
Donc, la série dLPNB1 est stationnaire.

3. Ajuster un modèle ARIMA pour la série LPNB1 c'est-à-dire un modèle ARMA pour dLPNB1
acf(dLPNB1); pacf(dLPNB1)

Series dLPNB1 Series dLPNB1

1.0

0.2
0.8

0.1
0.6

Partial ACF
0.4
ACF

0.0
0.2

-0.1
0.0

-0.2
-0.2

0 5 10 15 5 10 15

Lag Lag

Selon les fonctions acf et pacf de la série différenciée dLPNB, nous pouvons proposer un modèle
AR(1) pour cette série.
Fit=armaFit(~arma(1,0),data=dLPNB1);summary(Fit)
Title: ARIMA Modelling ;
Call: armaFit(formula = ~arma(1, 0), data = dLPNB1)
Coefficient(s): Estimate Std. Error t value Pr(>|t|)
ar1 0.250415 0.114746 2.182 0.0291 *
intercept 0.006872 0.001661 4.138 3.51e-05 ***
sigma^2 estimated as: 0.0001105; log likelihood: 222.65;AIC Criterion: -439.29
Autre méthode
library(forecast); auto.arima(LPNB1)
Series: LPNB1 ; ARIMA(1,1,0) with drift
Coefficients: ar1 drift
0.2504 0.0069
s.e. 0.1147 0.0017
sigma^2 estimated as 0.0001105: log likelihood=222.65

12
 AIC=-439.29 AICc=-438.93 BIC=-432.5
4. Graphe et acf des résidus. Tests sur les residus
res=residuals(Fit);ts.plot(res);acf(res);

Series res

1.0
0.02

0.8
0.01

0.6
0.00

0.4
ACF
res

-0.01

0.2
-0.02

0.0
-0.03

-0.2
0 10 20 30 40 50 60 70 0 5 10 15

Time Lag

*) Test d’absence d’autocorrélation des résidus

Box.test(res, lag=20,type="Box-Pierce")
Box-Pierce test : data: res
X-squared = 14.2983, df = 20, p-value = 0.8151
**) Test d’homoscédasticité des résidus
library(FinTS); ArchTest(res)
Chi-squared = 7.5098, df = 12, p-value = 0.8222
p-value > 0.05→ hypothèse d’homoscédasticité acceptée
***) Test de normalité des résidus
Jarque.bera.test(res)
Jarque Bera Test:data: res
X-squared = 2.293, df = 2, p-value = 0.3177
p-value > 0.05→ normalité des résidus.
Nous concluons que les résidus res constituent un bruit blanc gaussien.
5. Prévisions pour l’année 1988
Fit1=arima(LPNB1,order=c(1,1,0));
Prev=predict(Fit1,4)
Année: Trimestre Prev$pred Prev$se
1988 : Q1 8.2814 0.0113
1988 : Q2 8.2846 0.0201
1988 : Q3 8.2861 0.0277
1988 : Q4 8.2867 0.0343

13
Annexe

Test de stationnarité de la série LPNB1 par le test de Dickey-Fuller

x1=ur.df(LPNB,lag=1,type='trend'); summary(x1)
# Augmented Dickey-Fuller Test Unit Root Test # Test regression trend
Value of test-statistic is: -2.1636 10.8195 2.4789
Critical values for test statistics: 1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
tc =-2.16> tau3=-3.99 →accepter H0 : rau=0 dans (1) présence de racine unitaire
Fc =2.47< phi3=8.43→accepter H0 : rau=0 beta=0 dans (1)
On passe au modèle (2)
x2=ur.df(LPNB,lag=1,type='drift'); summary(x2)
# Augmented Dickey-Fuller Test Unit Root Test # Test regression drift
Value of test-statistic is: -0.7932 13.7257
Critical values for test statistics: 1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
tc =-0.79 > tau2 =-3.46 →accepter H0 : rau=0 dans (1) présence de racine unitaire
Fc =13.72> phi1=6.52→refuser H0 : rau=0 alpha=0 dans (2)→ LPNB1est un I(1)+T

Exemple empirique d’un modèle SARIMA

Nous étudions la série chronologique du nombre de passagers (en milliers) par mois, qu’on
note NP, dans les transports aériens, de 1949 à 1960.
Cette série est disponible sous R (AirPassengers).
1. Représenter graphiquement la série NP . Ce processus présente-t-il une stabilité de la variance, une
tendance et une saisonnalité. Que peut-on dire sur la stationnarité de ce processus ?
2. Tracer le graphe de la série LNP=log(NP). Donner une analyse visuelle de cette série.
3. Appliquer la méthode des différences pour enlever la saisonnalité et la tendance. La série obtenue
ddLNP semble-t-elle stationnaire ? Utiliser le test de Dickey –Fuller pour confirmer votre réponse.
4. Tracer graphiquement les fonctions d’auto-corrélation simple et partielle de la série ddLNP. Donner
des commentaires. Identifier un modèle ARMA pour cette série.
5. Estimer les paramètres du modèle choisi.
6. Utiliser les tests de Box-Tierce , Arch et de Jarque-Bera pour la série des résidus obtenus.
7. Calculer les prévisions pour l’année 1961.

14
1. NP=AirPassengers ;

Le graphique de la série et celui de la fonction d’autocorrélation sont donnés par ,

ts.plot(NP,xlab="Temps",ylab="Nombre de voyageurs") ;
600
500
Nombre de voyageurs

400
300
200
100

1950 1952 1954 1956 1958 1960

Temps

acf(NP,40,"correlation")

Series NP
1.0
0.8
0.6
ACF

0.4
0.2
0.0
-0.2

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Lag

Graphiquement nous notons que la série présente une tendance, une saisonnalité et une augmentation
de la variance. La fonction d’autocorrélation ne décroit pas rapidement. Ce qui signifie que la série est
non stationnaire.
2. Pour la stabilité de la variance nous considérons le logarithme népérien de la série brute. le graphe
de la serie LNP=log(NP) est donné par,
LNP=log(NP) ;ts.plot(LNP,xlab="Temps",ylab="LNP") ;

15
 6.5
6.0
LNP

5.5
5.0

1950 1952 1954 1956 1958 1960

Temps

Nous notons une stabilité de la variance de la série LNP mais elle présente toujours une tendance et
une saisonnalité.
3. Pour l’élimination de la saisonnalité de la série LNP nous considérons la série,
dLNP=diff(LNP,12)
Le graphe de cette série ainsi celui de la fonction fac sont donnés par,
ts.plot(dLNP,xlab="Temps",ylab="dLNP") ; acf(dLNP,40,"correlation")
Series dLNP
1.0
0.3

0.8
0.6
0.2

0.4
dLNP

ACF
0.1

0.2
0.0
0.0

-0.2

1950 1952 1954 1956 1958 1960 0.0 0.5 1.0 1.5 2.0 2.5 3.0

Temps Lag

Nous notons graphiquement que la série désaisonnalisée dLNP est non stationnaire. Le test de
stationnarité de la série dLNP est abordé en annexe.
Pour l’élimination de la tendance de la série dLNP nous considérons la série, ddLNP=diff(dLNP).
Le graphe de cette série est obtenu par, ts.plot(ddLNP,xlab="Temps",ylab="ddLNP") ;

16
 0.15
0.10
0.05
ddLNP

0.00
-0.05
-0.10
-0.15

1950 1952 1954 1956 1958 1960

Temps

La série ddLNP parait stationnaire. Vérifions la stationnarité de la série ddLNP par le test de Dickey-
Fuller. avec un retard p=11 qu’on détermine à partir de la fonction d’autocorrélation partielle .
pacf(ddLNP, 40)

Series ddLNP
0.6
0.4
0.2
Partial ACF

0.0
-0.2
-0.4

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Lag

x1=ur.df(ddLNP,lags=11,type='trend'); summary(x1);
# Augmented Dickey-Fuller Test Unit Root Test # Test regression trend
Estimate Std. Error t value Pr(>|t|)
tt -0.0001806 0.0001011 -1.785 0.076784 .
Value of test-statistic is: -4.9398 8.2389 12.353
Critical values for test statistics: 1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
- La valeur ( t c  t ˆ  4.93 ) est inférieure à t tab  tau 3  3.99 →rejeter H0 : ρ = 0 dans le modèle (1)

- La p-value du paramètre β (0.076) supérieur à 0.05→ accepter H 0 :β = 0 dans le modèle (1)

Maintenant, nous appliquons l’instruction :
x2=ur.df(ddLNP,lags=11,type='drift'); summary(x2)
pour obtenir les résultats suivants :
# Augmented Dickey-Fuller Test Unit Root Test # Test regression drift

17
(Intercept) 0.034084 0.008225 4.144 6.45e-05 ***
Value of test-statistic is: -4.5964 10.5686
Critical values for test statistics: 1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
La statistique ( t c  t ˆ  4.59 ) est inférieure à t tab  tau 2  3.46 →rejeter H0 : ρ = 0 dans le modèle (2).

- La p-value du paramètre α (6.45*10-5 ) inférieure à 0.05→ rejeter H0 : α = 0 dans le modèle (2) . Donc, la
série ddLNP est de la forme I(0)+c
Nous déduisons la stationnarité de la série ddLNP.
4. Les fonctions d’autocorrélation simple et partielle de la série ddLNP sont données par,
acf(ddLNP,40,"correlation") ; pacf(ddLNP,40,"correlation")

Series ddLNP Series ddLNP

1.0

0.2
0.8

0.1
0.6

0.0
Partial ACF
0.4
ACF

-0.1
0.2
0.0

-0.2
-0.2

-0.3
-0.4

0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0

Lag Lag

Nous notons que, les coefficients ACF sont significatifs aux retards 1, 3 (?) et 12 , 23(?). et que les
coefficients PACF sont significatifs aux retards 1, 3(?), 9, 12
Nous suggérons à la série brute LN un modèle SARIMA(p,d,q)x(P,D,Q)12 avec d=1, D=1,
pmax=2,qmax=2,Pmax=1 et Qmax=1.
Nous utilisons la fonction auto.arima( ) du package forecast pour le meilleur modèle
library(forecast)
best.mod=auto.arima(LNP,d=1,D=1,max.p=2,max.q=2,max.P=1,max.Q=1)
best.mod; Series: LNP ; ARIMA(0,1,1)(0,1,1)[12]
5. Estimation des paramètres du modèle choisi
Fit=arima(LNP,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
Series: LNP ;ARIMA(0,1,1)(0,1,1)[12]
Coefficients: ma1 sma1
-0.4018 -0.5569
s.e :. 0.0896 0.0731
sigma^2 estimated as 0.001348: log likelihood=244.7
AIC=-483.4 AICc=-483.21 BIC=-474.77

18
6. Graphe et acf des résidus. Tests sur les résudus
Les résidus sont donnés par , res=residuals(Fit).
A partir des instructions, ts.plot(res); acf(res,40,"correlation"), nous obtenons
le graphe et la fonction acf des résidus.
Series res

1.0
0.10

0.8
0.05

0.6
0.00

ACF
res

0.4
-0.05

0.2
0.0
-0.10

-0.2

1950 1952 1954 1956 1958 1960 0.0 0.5 1.0 1.5 2.0 2.5 3.0

Time Lag

Test d’absence d’autocorrélation des résidus

Box.test(res, lag=40,type="Box-Pierce")
X-squared = 34.7935, df = 40, p-value = 0.7033
Les résidus ne sont pas autocorréles car p-value=0.70 > 0.05
Test d’homoscédasticité des résidus
library(FinTS) ;ArchTest(res)
data: res ; Chi-squared = 14.859, df = 12, p-value = 0.2493
p-value =0.24 > 0.05→ hypothèse d’homoscédasticité acceptée.
Test de normalité des résidus
library(tseries)) ;jarque.bera.test(res)
X-squared = 5.2265, df = 2, p-value = 0.0733

Nous acceptons l’hypothèse de normalité car, JQ  5.22   2 (0.05)  5.99

(2)
ou p-value =0.07 > 0.05
7. Calculer les prévisions pour l’année 1961.
p=predict(Fit,n.ahead=12);
Valeurs prévues:
prev1961=exp(p$pred);
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
450.4 425.7 479.0 492.4 509.1 583.3 670.0 667.1 558.2 497.2 429.9 477.2
Intervalle de prevision
Les bornes inférieures et supérieures des intervalles de prévision sont données par,
A=exp(p$pred-1.96*p$se) ;B= exp(p$pred+1.96*p$se) ;

19
Jan Fev Mar Avr Mai Juin Juillet Aout Sep Oct Nov Dec
419.1 391.5 435.9 443.9 455.0 517.3 589.7 583 484.6 428.9 368.5 406.7
484.0 462.9 526.3 546.2 569.5 657.8 761.2 763.3 643 576.4 501.4 560
Graphe des prévisions de la série du nombre de passagers
NPprev=ts(c(NP,prev1961),start=1949,frequency=12);plot(NPprev,type="l",lty=1, col="blue")
lines(A, lty=2,col="red");lines(B, lty=2,col="red")
600
500
NPprev

400
300
200
100

1950 1952 1954 1956 1958 1960 1962

Time

Annexe
Test de stationnarité de Dickey-Fuller pour la série dLNP.
Utilisons le test de Dickey-Fuller avec un retard p=13 qu’on détermine à partir de la fonction d’autocorrélation
partielle.
x1=ur.df(dLNP,lags=13,type='trend'); summary(x1)
Puisque tc  t ˆ  2.87  ttab  tau3  3.99, nous acceptons H0 : ρ = 0 pour un seuil de 1% dans le modèle (1).

D’autre part, l’hypothèse (β, ρ) = (0, 0) est acceptée dans le modèle (1) car phi3c= 4.12 < phi3t = 8.43 au seuil de 1% .
Donc on accepte l’hypothèse que la série a une racine unité et pas de tendance au un seuil de 1%.
x2=ur.df(dLNP,lags=13,type='drift'); summary(x2)

L’hypothèse ρ = 0 est acceptée pour un seuil de 1% dans le modèle (2) car t c  t ˆ  2.393  t tab  tau 2  3.46.
D’une manière similaire, nous dirons que l’hypothèse (α, ρ) = (0, 0) est acceptée dans le modèle (2) car phi1c =6.52 < phi1t
=6.52 au seuil de 1% .Donc on accepte l’hypothèse que la série a une racine unité et ne possède pas de constante pour un
niveau de 1%.
x3=ur.df(dLNP,lags=13,type='none'); summary(x3);

L’hypothèse ρ = 0 est acceptée pour un seuil de 1% dans le modèle (3) car t c  t ˆ  0.92  t tab  tau1  2.58.
Nous pouvons conclure que pour un seuil de 1% la série désaisonnalisée dLNP présente une racine unitaire et ne possède pas
de tendance ni de constante. La série est intégrée d’ordre 1, c'est-à-dire I(1) et donc nous la rendons stationnaire en
considérant la série différenciée d’ordre 1. C'est-à-dire, la série ddLNP qui est égale à diff(dLNP).

20