Séries Temporelles 1 (2A) - ENSAI

TP 1 : Premiers pas et quelques expériences numériques

Octobre 2023

Exercice 1
On étudie la série temporelle nottem des températures mensuelles enregistrées à Nottingham Castle entre
janvier 1920 et décembre 1939. Cette série est disponible sous R :
[Link] <- nottem
summary([Link])

Min. 1st Qu. Median Mean 3rd Qu. Max.

31.30 41.55 47.35 49.04 57.00 66.50

1.- Regarder la structure de ces données, puis tracer la série à l’aide de la fonction [Link]. Observe-t-on
une tendance ? une saisonnalité ?
Voici la répresentation graphique en format de série temporelle ts :
[Link]([Link])
30 35 40 45 50 55 60 65
[Link]

1920 1925 1930 1935 1940

Time
# Tracez le même graphe avec plot([Link]). Remarquez-vous une différence ?

A simple vue, on n’observe pas de tendance (globale) notable. On pourrait valider cela avec un test. Par
exemple,
𝐻0 ∶ {𝑋𝑡 = 𝑎𝑡 + 𝑏 , 𝑎 ≠ 0} (existence d’une tendance linéaire)

1
vs
𝐻1 ∶ {𝑋𝑡 = 𝑏} (pas de tendance linéaire)

Pour cela, il suffit d’utiliser la commande lm et voir le résultat du test donnée par la commande summary.
mois <- 1:length([Link])
mod0 <- lm([Link] ~ mois)
summary(mod0)

Call:
lm(formula = [Link] ~ mois)

Residuals:
Min 1Q Median 3Q Max
-17.675 -7.575 -1.620 8.028 17.882

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 48.301984 1.111112 43.472 <2e-16 ***
mois 0.006121 0.007994 0.766 0.445
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.58 on 238 degrees of freedom

Multiple R-squared: 0.002458, Adjusted R-squared: -0.001734
F-statistic: 0.5864 on 1 and 238 DF, p-value: 0.4446
On ne peut pas rejeter 𝐻0 , on valide donc qu’on n’a pas de tendance statistiquement significative. Cependant,
en regardant les données, on observe bien une saisonalité (qui n’est pas étonnant, vue la nature des données)
Remarque : Vous pouvez vérifier cela en utilisant le periodogram et/ou un test de saisonalité (non vu en
cours)
TSA::periodogram([Link])
abline(v=1/12, col="red", lty=2)

2
 15000
10000
Periodogram

5000
0

0.0 0.1 0.2 0.3 0.4 0.5

Frequency
Ce graphe montre bien que nous avous une période de 𝑃 = 12, car 𝑓𝑟𝑒𝑞 = 1/𝑃

2.- On tente de la désaisonnaliser par régression sur les fonctions trigonométriques, 𝑡 → cos ( 2𝜋𝑗12 𝑡) pour
2𝜋𝑗
𝑗 = 1, 2, … , 6 et 𝑡 → sin ( 12 𝑡) pour 𝑗 = 1, … , 5. Pour cela, on crée une base trigonométrique à l’aide du
code suivant :
freq <- matrix(1:6, nrow=1)/12
time <- matrix(1:length([Link]), ncol = 1)
B <- cbind(cos(2*pi*time%*%freq), sin(2*pi*time%*%freq))[,-12]
B <- [Link](B)
colnames(B) <- c('cos1','cos2','cos3','cos4','cos5','cos6','sin1','sin2','sin3',
'sin4','sin5')
library(knitr)
knitr::kable(head(B), digits = 4) #Si vous voulez regarder les premières 6 lignes du tableau B

cos1 cos2 cos3 cos4 cos5 cos6 sin1 sin2 sin3 sin4 sin5
0.866 0.5 0 -0.5 -0.866 -1 0.500 0.866 1 0.866 0.500
0.500 -0.5 -1 -0.5 0.500 1 0.866 0.866 0 -0.866 -0.866
0.000 -1.0 0 1.0 0.000 -1 1.000 0.000 -1 0.000 1.000
-0.500 -0.5 1 -0.5 -0.500 1 0.866 -0.866 0 0.866 -0.866
-0.866 0.5 0 -0.5 0.866 -1 0.500 -0.866 1 -0.866 0.500
-1.000 1.0 -1 1.0 -1.000 1 0.000 0.000 0 0.000 0.000

3.- Effectuer la régression sur la base créée en utilisant la fonction lm puis faire un résumé (summary) des
résultats obtenus.
L’idée est de réaliser ici une régression (parfois appelé “régression harmonique”) :
6 5
2𝜋𝑗 2𝜋𝑗
𝑋𝑡 = 𝜇 + ∑ 𝛼𝑗 cos ( 𝑡) + ∑ 𝛽𝑗 sin ( 𝑡) + 𝜖𝑡
𝑗=1
12 𝑗=1
12

3
afin d’obtenir une approximation de la partie saisonnière 𝑆𝑡 de la decomposition additive de la série (𝑋𝑡 ),
i.e., 𝑋𝑡 = 𝑆𝑡 + 𝑇𝑡 + 𝜖𝑡 .
data_base <- [Link](cbind(B,[Link]))
colnames(data_base)[12] <- 'temperature'
mod0 <- lm(temperature~. , data=data_base)
summary(mod0)

Call:
lm(formula = temperature ~ ., data = data_base)

Residuals:
Min 1Q Median 3Q Max
-7.890 -1.369 0.285 1.405 6.270

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.03958 0.14942 328.207 < 2e-16 ***
cos1 -9.24092 0.21131 -43.732 < 2e-16 ***
cos2 -0.08083 0.21131 -0.383 0.7024
cos3 -0.06417 0.21131 -0.304 0.7617
cos4 0.02167 0.21131 0.103 0.9184
cos5 0.05009 0.21131 0.237 0.8128
cos6 -0.19542 0.14942 -1.308 0.1922
sin1 -6.94091 0.21131 -32.848 < 2e-16 ***
sin2 1.49822 0.21131 7.090 1.66e-11 ***
sin3 0.34333 0.21131 1.625 0.1056
sin4 0.36517 0.21131 1.728 0.0853 .
sin5 0.14174 0.21131 0.671 0.5030
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.315 on 228 degrees of freedom

Multiple R-squared: 0.9304, Adjusted R-squared: 0.9271
F-statistic: 277.3 on 11 and 228 DF, p-value: < 2.2e-16
On observe que les éléments significatifs sont cos1 (𝛼1 ), sin1 (𝛽1 ), sin2 (𝛽2 ) et l’intercept (𝜇).

4.- Recommencer en ne gardant que les éléments significatifs de la base.

mod1 <- lm(temperature~cos1+sin1+sin2, data=data_base)
summary(mod1)

Call:
lm(formula = temperature ~ cos1 + sin1 + sin2, data = data_base)

Residuals:
Min 1Q Median 3Q Max
-8.4056 -1.4098 0.3104 1.3961 6.0013

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.0396 0.1494 328.143 < 2e-16 ***

4
cos1 -9.2409 0.2113 -43.724 < 2e-16 ***
sin1 -6.9409 0.2113 -32.841 < 2e-16 ***
sin2 1.4982 0.2113 7.089 1.56e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.315 on 236 degrees of freedom

Multiple R-squared: 0.928, Adjusted R-squared: 0.9271
F-statistic: 1014 on 3 and 236 DF, p-value: < 2.2e-16

5.- Créer une série temporelle contenant la partie saisonnière obtenue sur cette base puis une contenant les
résidus de la régression. On fera appel aux commandes $fitted et $residuals.
• Voici la série de saisonalité 𝑆𝑡 :
#On utilise la commande $fitted :
S_t <- mod1$[Link]
s_t <- ts(S_t, frequency = 12, start = c(1920,1)) #Pour obtenir la série placée
#sur les dates des observations [Link]

[Link]([Link]) #La série originale

lines(s_t, col="red", lty=2) #La série S_t est la rouge pointillée
30 35 40 45 50 55 60 65
[Link]

1920 1925 1930 1935 1940

Time
• Voici la série des résidus 𝜖𝑡 :
E_t <- mod1$residuals
e_t <- ts(E_t, frequency = 12, start = c(1920,1)) #Pour obtenir la série placée
#sur les dates des observations [Link]

plot(e_t)

5
 5
0
e_t

−5

1920 1925 1930 1935 1940

Time
Bonus : vous pouvez vérifier si ces résidus estimés sont par exemple centrés, gaussiens, etc. Pour cela, vous
pouvez tracer des histogrammes, qq-plots, etc.

Vous pourriez également trouver la partie saisonnière et les résidus d’une manière “manuelle”, c’est-à-dire,
sans utiliser les commandes $[Link] et $residuals :
## Solution manuelle

#partie saisonniere :
S_t2 <- matrix(c(rep(1, length(time)),cos(pi*time/6), sin(pi*time/6), sin(pi*time/3)),
ncol=4, byrow = FALSE)%*%[Link](mod1$coefficients)
s_t2 <- ts(S_t2, frequency = 12, start = c(1920,1))
[Link]([Link])
lines(s_t2, col="red", lty=2)

6
 30 35 40 45 50 55 60 65
[Link]

1920 1925 1930 1935 1940

Time
#serie des residus :
E_t2 <- [Link] - S_t2
e_t2 <- ts(E_t2, frequency = 12, start = c(1920,1))
[Link](e_t2)
5
Series 1

0
−5

1920 1925 1930 1935 1940

Time

6.- Tracer les trois séries dans une même fenêtre.

par(mfrow=c(3,1))
[Link]([Link], main="Données observées", ylab = "observed")
[Link](s_t, main="Partie saisonière", ylab="seasonal")

7
[Link](e_t, main="Résidus", ylab="random")

Données observées
observed

30 60

1920 1925 1930 1935 1940

Time

Partie saisonière
seasonal

40 60

1920 1925 1930 1935 1940

Time

Résidus
random

−5

1920 1925 1930 1935 1940

Time

7.- Filtrer la série de départ pour éliminer la saisonnalité en utilisant le filtre de différenciation saisonnière
via la fonction diff. Superposer la série ainsi obtenue avec les résidus précédents. Commenter.
On a vu en cours que le filtre Δ𝑝 = (𝐼 − 𝐵𝑝 ) absorbe les saisonnalités de période 𝑝. Sachant que notre série
a une composante saisonnière de période 𝑝 = 12 on applique la commande suitante :
D_12 <- diff([Link], lag=12) #lag=p (la période de la série)

# Essayez de tester plusieurs lags pour observer si la série résultante conserve

# ou non une saisonnalité.

par(mfrow=c(1,1))
plot(D_12, col="red", lty=2, ylab="residuals")
lines(e_t)

8
 10
5
residuals

0
−5
−10

1925 1930 1935 1940

Time
Sur le graphe, la trajectoire rouge correspond à la série filtrée et la trajectoire noire à la série des résidus
estimés précedement. On peut constater une variance plus élévée dans le cas de la série filtrée et c’est dû à
l’effet du filtre. En effet, si

𝑋𝑡 = 𝑚 + 𝑆𝑡 + 𝜖𝑡 ⟹ 𝜖′𝑡 ∶= Δ12 (𝑋𝑡 ) = 𝜖𝑡 − 𝜖12 .

On peut observer cette dispersion en regardant les densités (empiriques) :

#hist(e_t, probability = TRUE, breaks = 12)
plot(density(e_t), col="black", main = "densités de e_t et e'_t")
lines(density(D_12), col="red")

9
 densités de e_t et e'_t
0.20
0.15
Density

0.10
0.05
0.00

−10 −5 0 5

N = 240 Bandwidth = 0.6298

8.- Même question en utilisant cette fois la moyenne mobile arithmétique d’ordre 13 modifiée 𝑀6∗ (𝐵). On
utilisera la fonction filter qui prend en argument la série à filter et un vecteur contenant les coefficients 𝜃𝑖
de la moyenne mobile à appliquer.
∗ 1 𝑚−1 1
On sait que 𝑀𝑚 (𝐵) = ∑ 𝐵−𝑘 + (𝐵−𝑚 + 𝐵𝑚 ) est un filtre qui absorbe la saisonnalité de
2𝑚 𝑘=−(𝑚−1) 4𝑚 ∗
période 𝑝 = 2𝑚. Ici, 𝑝 = 12. Les coefficients de 𝑀6 (𝐵) sont donc dans le vecteur de dimension 13 :

1 1 1 1
(𝜃−6 , 𝜃−5 , … , 𝜃5 , 𝜃6 ) = ( , ,…, , )
24 12 12 24

M_6_modif <- filter([Link],

filter = c(1/24, rep(1/12, 11), 1/24) #les coefficients du filtre
)
#plot(M_6_modif, col="red", lty=1)
[Link]([Link])
lines(M_6_modif, col="red", lty=1)

10
 30 35 40 45 50 55 60 65
[Link]

1920 1925 1930 1935 1940

Time
On voit bien qu’il n’y a plus de saisonnalité. À différence du filtre Δ12 , ici, la série filtré est centrée autour
de la moyenne globale. (Exercice, Pourquoi ?)

9.- Pour terminer, utiliser la fonction decompose. Que fait-elle? Faire un graphe de sa sortie. Commenter.
decomposition <- decompose([Link], type="additive")
plot(decomposition)

Decomposition of additive time series

observed
60
45
30
51
trend
49
10 47
random seasonal
0
−10
2
−2
−6

1920 1925 1930 1935 1940

Time

11
Regardez le suivant :
sum([Link](decomposition$trend - M_6_modif))

[1] 0
Cela montre que la série “trend” de decompose est 𝑀6∗ (𝐵)𝑋𝑡 .
• Comment sont-elles calculés les séries seasonal et random (ou residuals) dans la fonction decompose?
(Discuté en TP)

Exercice 2
Soient (𝑋𝑡 )𝑡 et (𝑌𝑡 )𝑡 deux processus MA définis pour tout 𝑡 ∈ ℤ par

1 1
𝑋𝑡 = 3 + 𝜖𝑡 + 2𝜖𝑡−1 𝑒𝑡 𝑌𝑡 = 𝜂𝑡 − 𝜂𝑡−1 + 𝜂𝑡−2
2 3
où (𝜖𝑡 )𝑡 et (𝜂𝑡 )𝑡 sont deux bruits blancs indépendants de variance respective 𝜎𝜖2 = 1 et 𝜎𝜂2 = 2. On pose
𝑍𝑡 = 𝑋𝑡 + 𝑌𝑡 .
1.- Simuler une trajectoire de 𝑋, de 𝑌 et de 𝑍, de taille 𝑛 = 1000.
[Link](666)
n <- 1000
m <- n + 2
epsilon <- rnorm(m, sd=1)
eta <- rnorm(m, sd=sqrt(2))
X <- 3 + epsilon[-1] + 2*epsilon[-m]
X <- X[1:n] #MA(1)
Y <- eta[-c(1,2)]-1/2*eta[-c(1,m)]+1/3*eta[-c(m-1,m)]
Y <- Y[1:n] #MA(2)
Z <- X + Y

par(mfrow=c(2,2))
[Link](X)
[Link](Y)
[Link](Z)

12
 8

4
X

−4 0
2
−4

0 200 400 600 800 1000 0 200 400 600 800 1000

Time Time
5
Z

−5

0 200 400 600 800 1000

Time
2.- A l’aide des fonctions acf et pacf, tracer les fonctions d’autocorrélation et d’autocorrélation partielle
empirique de chacun des trois processus.
par(mfrow=c(1,2))
pacf(X)
acf(X)

Series X Series X
1.0
0.3

0.8
0.2

0.6
Partial ACF

ACF
0.1

0.4
0.0

0.2
0.0
−0.2

0 5 10 20 30 0 5 10 20 30

Lag Lag

13
par(mfrow=c(1,2))
pacf(Y)
acf(Y)

Series Y Series Y

1.0
0.1

0.5
Partial ACF

−0.1

ACF

0.0
−0.3

−0.5
−0.5

0 5 10 20 30 0 5 10 20 30

Lag Lag
par(mfrow=c(1,2))
pacf(Z)
acf(Z)

14
 Series Z Series Z

1.0
0.06

0.8
0.02

0.6
Partial ACF

ACF

0.4
−0.02

0.2
−0.06

0.0
0 5 10 20 30 0 5 10 20 30

Lag Lag

Exercice 3
Soit le processus AR(1) défini, pour tout 𝑡 ∈ ℤ, par

𝑋𝑡 = 𝜇 + 𝜙𝑋𝑡−1 + 𝜖𝑡

où (𝜖𝑡 )𝑡 est un brut blanc de variance 𝜎2 .

1.- Simuler une trajectoire de ce processus de taille 𝑛 = 1000 pour 𝜇 = 1, 𝜎2 = 0.5 et 𝜙 = 0.45, puis pour
𝜙 = 0.95 et enfin pour 𝜙 = 1.
Pour que cela soit automatique, on construit d’abord la function suivante :
AR1 <- function(mu,phi,epsilon,size){
X <-c(epsilon[1]) + mu
for (t in 2:length(epsilon)){
X[t] <- mu + phi*X[t-1] + epsilon[t]
}
return(X[1:size])
}

et on utilise cette fonction pour simuler chaque série temporelle :

[Link](1001)
epsilon1 <- rnorm(1010, sd=sqrt(0.5))
X1 <- AR1(mu=1, phi=0.45, epsilon = epsilon1, size=1000)
X2 <- AR1(mu=1, phi=0.95, epsilon = epsilon1, size=1000)
X3 <- AR1(mu=1, phi=1, epsilon = epsilon1, size=1000)

2.- Pour chaque valeur de 𝜙 représenter graphiquement la trajectoire simulée et commentez-les.

15
par(mfrow=c(2,2))
[Link](X1)
[Link](X2)
[Link](X3)

15 25
4
X1

X2
2

5
0

0 200 400 600 800 1000 0 200 400 600 800 1000

Time Time
600
X3

0 200 400 600 800 1000

Time
La première série semble être stationnaire. La troisième série est clairement non-stationnaire. La deuxième
série semble être entre les deux.
Remarque : Dans la théorie, la deuxième série est encore stationnaire. En effet, la série

𝑋𝑡 = 𝜇 + 𝜙𝑋𝑡−1 + 𝜖𝑡

est stationnaire si |𝜙| < 1. La valeur 𝜙 = 0.95 est presque à la frontière de la non-stationnarité.
3.- A l’aide de la fonction acf, tracer la fonction d’autocorrelation empirique du processus. Dans quels cas
est-il stationnaire ?
par(mfrow=c(1,3))
acf(X1)
acf(X2)
acf(X3)

16
 Series X1 Series X2 Series X3
1.0

1.0

1.0
0.8

0.8

0.8
0.6

0.6

0.6
ACF

ACF

ACF
0.4

0.4

0.4
0.2

0.2

0.2
0.0

0.0

0.0
0 5 10 20 30 0 5 10 20 30 0 5 10 20 30

Lag Lag Ces Lag

graphiques confirment ce qu’on a dit prècedement. C’est-à-dire, que les séries 1 et 2 sont stationnaires
et la 3eme non-stationnaire.
4.- Dans les cas stationnaires, estimer pour chaque 𝑘 ∈ {10, … , 𝑛}, 𝜇, 𝜙 et 𝜎2 à partir de la trajectoire
𝑥1 , … , 𝑥𝑘 (on utilisera la moyenne empirique et les équations de Yule-walker). Comparer graphiquement les
estimations avec les vraies valeurs des paramètres en fonction de la valeur de 𝑘.
Ici, les équations de Yule-Walker sont en effet une seule équation :

𝜙 = 𝜌(0)−1 𝜌(1) = 𝜌(1)

De cette manière, nous pouvons étudier le comporement des estimateurs de 𝜇, 𝜙 et 𝜎2 (𝜇𝑘̂ , 𝜙𝑘̂ , 𝜎̂𝑘 ) via la
fonction suivante :
[Link].AR1 <- function(X){
n <- length(X)
[Link] <-c(NA)
[Link] <- c(NA)
hat.sigma2 <- c(NA)
for (k in 10:n){
aux <- acf(X[1:k], [Link] = 2, plot = FALSE)
[Link][k] <- aux$acf[2] #Si l'on simplifie les equations de YW pour AR(1)
[Link][k] <- ([Link][k])*mean(X[1:k])
hat.sigma2[k] <- ([Link][k]^2)*var(X[1:k])
}
Parameters <- [Link](cbind([Link], [Link], hat.sigma2))
names(Parameters) <- c("phi", "mu", "sigma2")
return(Parameters)
}

• Pour la première série :

17
E1 <- [Link].AR1(X1)
par(mfrow=c(3,1))
plot(E1$phi, type="l", xlab="k")
abline(h=0.45, col="red")

plot(E1$mu, type="l")
abline(h=1, col="red")

plot(E1$sigma2, type = "l")

abline(h=0.5, col="red")
E1$phi

0.25

0 200 400 600 800 1000

k
E1$mu

1.0

0 200 400 600 800 1000

Index
E1$sigma2

0.50

0 200 400 600 800 1000

Index

• Pour la deuxième série :

E2 <- [Link].AR1(X2)
par(mfrow=c(3,1))
plot(E2$phi, type="l", xlab="k")
abline(h=0.95, col="red")

plot(E2$mu, type="l", xlab="k")

abline(h=1, col="red")

plot(E2$sigma2, type = "l", ylim = c(0, max([Link](E2$sigma2))))

abline(h=0.5, col="red")

18
E2$phi

0.60
0 200 400 600 800 1000

k
0.8 2.0
E2$mu

0 200 400 600 800 1000

k
E2$sigma2

0 3

0 200 400 600 800 1000

Index

5.- Pour chaque 𝑘 ∈ {10, … , 𝑛−1}, proposer un estimateur de la prévision de 𝑋𝑘+1 à partir de l’observation de
la trajectoire 𝑥1 , … , 𝑥𝑘 . Superposer la vraie trajectoire et sa prévision en fonction de 𝑘. Que constatez-vous
selon la valeur de 𝜙 ?
Prevision.AR1_1 <- function(X){
previ <- c(NA)
n <- length(X)
for (k in 10:n){
aux <- acf(X[1:k], [Link] = 2, plot = FALSE)
[Link] <- aux$acf[2]
[Link] <- ([Link])*mean(X[1:k])
previ[k+1] <- [Link] + [Link]*X[k] #Voir la fin du Chapitre 3
}
return(previ)
}

• Prévisions 𝑋̂ 𝑘 (1) pour la première série :

Prev.X1 <- Prevision.AR1_1(X1)
[Link](X1)
lines(Prev.X1, col="red", lty=1)

19
 4
3
X1

2
1
0

0 200 400 600 800 1000

Time
• Prévisions 𝑋̂ 𝑘 (1) pour la deuxième série :
Prev.X2 <- Prevision.AR1_1(X2)
[Link](X2)
lines(Prev.X2, col="red")
25
20
15
X2

10
5

0 200 400 600 800 1000

Time

Exercice 4
1 1
Soit (𝑋𝑡 )𝑡 un processus MA(2) de paramètres 𝑚 = 3, 𝜃1 = et 𝜃2 = − .
2 3
1.- Simuler une trajectoire pour 𝑡 = 1, … , 300.

20
[Link](134)
mu = 3
theta1 = 1/2
theta2 = -1/3
n = 300
bruit <- rnorm(n+2)
N = length(bruit)
X <- mu + bruit[-c(1,2)] + theta1*bruit[-c(1,N)] + theta2*bruit[-c(N-1,N)]
X <- X[1:n]
[Link](X)
5
4
3
X

2
1
0

0 50 100 150 200 250 300

Time
2.- Utiliser la fonction arima pour estimer les paramètres 𝑚, 𝜃1 et 𝜃2 par maximum de vraisemblance à partir
des 200 premières observations.
ma2 <- arima(X[1:200], order = c(0,0,2))
ma2$coef

ma1 ma2 intercept

0.5493713 -0.2881389 2.9393657
3.- Ecrire un algorithme pour prédire les valeurs de 𝑋𝑡 pour 𝑡 = 200, 201, … , 300 et les comparer graphique-
ment avec les vraies valeurs observées.
̄
Il faut calculer d’abord 𝑃𝐹𝑇 (𝑋𝑇 +ℎ ), où 𝑇 = 200 et ℎ = 1, 2, ..., 100. 𝐹𝑇 = Vect(𝑋 𝑇 , 𝑋𝑇 +1 , …). Calcul
fait (ou à faire) en TP. Le reste est facile à coder
Sinon, on peut réaliser le calcul avec la librairie forecast :
library(forecast)
prevision_h <- forecast(ma2, h=100)
plot(prevision_h)
lines(ts(X[201:300], start = 201, frequency=1), col="red")

21
5
4
3
2
1
0 Forecasts from ARIMA(0,0,2) with non−zero mean

0 50 100 150 200 250 300

Pas très intéressant comme prévision. Exercice : bricolez une solution plus réaliste !

22