Importation des données

Data_poux = read.csv2("Data_poux_mensuel_log_2010_2024.txt") #chargement des donnees
Data_poux$Date=as.factor(Data_poux$Date)

Statistiques

Dimension du tableau

dim(Data_poux)
## [1] 174   2

Résumé statistique de la variable étudiée

summary(Data_poux$Indice)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.800   4.510   4.640   4.656   4.795   5.410

Dates

Data_poux$Date = as.Date(Data_poux$Date, format = "%d:%m:%Y")

Data_poux$mois = format(x=Data_poux$Date, format = "%m")
Data_poux$annee = format(Data_poux$Date, format = "%Y")

Tableau de Buys-Ballot

library(tidyr)
## Warning: le package 'tidyr' a été compilé avec la version R 4.4.2
Data_poux_buys_ballot = pivot_wider(Data_poux, id_cols = annee, names_from = mois, values_from = Indice)

Série temporelle

Data_poux_ts = ts(Data_poux$Indice, start=c(2010,1), end=2024, freq=12)
# Zoom sur la période de janvier 2010 à décembre 2019
zoom_Data_poux_ts = window(x=Data_poux_ts,
                        start = c(2010, 1),
                        end = c(2019, 12),
                        freq=12)

# Reste de la période
zoom_Data_poux_ts_apres = window(x=Data_poux_ts,
                        start = c(2020, 1),
                        freq=12)

Chronogramme

#graphique
plot(zoom_Data_poux_ts,
     xlab = "Années",
     ylab = "Indice",
     main = "Amélie Condette",
     xlim=c(2010,2024),
     ylim=c(3.8,5.4))

points(zoom_Data_poux_ts, 
       pch=20,
       col="darkblue")
lines(zoom_Data_poux_ts_apres, col="grey", type='o', pch=20)

Sur cet autocorrélogramme, on peut observer une saisonnalité ainsi qu’une tendance linéaire.

Autocorrélogramme

acf(zoom_Data_poux_ts, lag.max=60, main="Amelie Condette")

Cet autocorrélogramme montre une saisonnalité de période 12. ## Visualisation par année (courbes empilées)

matplot(t(Data_poux_buys_ballot[,-1]), # transposé :  1 courbe <=> 1 année
        type='b', #ligne avec des points
        xlab="Mois", #Nom axe des abscisses 
        ylab="Indice", #Nom axe des ordonnées
        col=topo.colors(15), # une couleur par année
        pch=20, #forme des points
        lty=1) #type de ligne

Ces courbes empilées confirment la saisonnalité pour chaque année avec un pic au mois d’août.

Décomposition de la série

#Décomposition de la série
decomposition_zoom_Data_poux_ts = decompose(x=zoom_Data_poux_ts,
          type = "additive")

# Visualisation de la décomposition
plot(decomposition_zoom_Data_poux_ts,
     col = "darkgreen")

D’après les graphiques précédents, nous avons un modèle additif de type \(x_t=s_t\).

Etude de la tendance et/ou saisonnalité

Méthode paramétrique

Estimation des composantes

# Tendance estimée
T= length(zoom_Data_poux_ts)
t=1:T
ajust_tend = lm(zoom_Data_poux_ts~t)
ajust_tend_ts = ts(ajust_tend$fitted.values,  start=2010, freq=12)

plot(zoom_Data_poux_ts, xlab="Années", ylab="Indice", col="black", main="Tendance estimée")
lines(ajust_tend_ts, col="blue", pch=20)

# Retranchement de la série
Data_poux_sans_tend = zoom_Data_poux_ts - ajust_tend$fitted.values
plot(Data_poux_sans_tend, xlab="Années", ylab="Indice", col="black", main="Série sans tendance")

# Estimation de la saisonnalité
S=12
sint = sin(2*pi*t/S)
cost = cos(2*pi*t/S)

ajust_sais = lm(Data_poux_sans_tend~sint+cost)
ajust_sais_ts = ts(ajust_sais$fitted.values,  start=2010, freq=12)

plot(Data_poux_sans_tend, xlab="Années", ylab="Indice", col="black", main="Saisonnalité estimée")
lines(ajust_sais_ts, col="blue", pch=20)

Prévisions pour des valeurs futures de la série (horizon précis)

# Prédiction de la tendance
N=length(zoom_Data_poux_ts_apres)
n=(T+1):(T+N)
nd_tend=data.frame(t=n)
pred_tend=predict(ajust_tend, newdata=nd_tend)
# Prédiction de la saisonnalité
N=length(zoom_Data_poux_ts_apres)
n=(T+1):(T+N)
nd_sais=data.frame(sint=sin(2*pi*n/S),
              cost=cos(2*pi*n/S)
              )
pred_sais=predict(ajust_sais, newdata=nd_sais)
# Prédiction de la série
predserie=ts(pred_tend+pred_sais,
             start=c(2020,1),
             freq=12)
estim_serie = ts(ajust_tend$fitted.values+ajust_sais$fitted.values, start=2010, freq=12)

plot(zoom_Data_poux_ts, xlab="Années", ylab="Indice", col="black", 
     xlim=c(2010, 2024), type="o", pch=20, main="Amelie Condette", ylim=c(3.8,5.4))
lines(zoom_Data_poux_ts_apres, col="grey", type="o", pch=20)
lines(estim_serie, col="blue", pch=15)
lines(predserie, col="red", type="o", pch="*")

Méthode non paramétrique

Estimation de la tendance à l’aide de moyennes mobiles

# Définition de l'ordre
p=12

# Estimation
Data_poux_avant_tend = filter(zoom_Data_poux_ts, c(1/(2*p), rep(1/p, p-1), 1/(2*p)))
Data_poux_avant_tend
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2010       NA       NA       NA       NA       NA       NA 4.795000 4.795000
## 2011 4.782083 4.771667 4.757083 4.747917 4.740833 4.737917 4.738750 4.730000
## 2012 4.718750 4.721250 4.726250 4.733750 4.742917 4.749583 4.752917 4.761667
## 2013 4.782917 4.786250 4.793750 4.802917 4.812917 4.823333 4.839583 4.859167
## 2014 4.828333 4.817500 4.800000 4.779583 4.759583 4.740000 4.715833 4.685417
## 2015 4.695833 4.698333 4.699167 4.696250 4.694583 4.700417 4.708750 4.717917
## 2016 4.681250 4.672083 4.665417 4.662917 4.660417 4.653750 4.642500 4.637500
## 2017 4.701667 4.718333 4.726250 4.729583 4.730417 4.726667 4.729583 4.728333
## 2018 4.615417 4.587500 4.571250 4.559167 4.546250 4.539167 4.527500 4.512500
## 2019 4.579583 4.595833 4.608750 4.620417 4.633750 4.643750       NA       NA
##           Sep      Oct      Nov      Dec
## 2010 4.794583 4.791250 4.787917 4.785417
## 2011 4.718750 4.717917 4.718333 4.717917
## 2012 4.771250 4.772500 4.776250 4.782500
## 2013 4.867917 4.864167 4.855417 4.840417
## 2014 4.670833 4.675417 4.682500 4.690833
## 2015 4.717917 4.707917 4.696250 4.686667
## 2016 4.647500 4.662917 4.675417 4.686667
## 2017 4.708333 4.682500 4.660000 4.640833
## 2018 4.516667 4.535417 4.555000 4.568333
## 2019       NA       NA       NA       NA
# Graphique
plot(zoom_Data_poux_ts, 
     xlab = "Années",
     ylab = "Indice")
lines(Data_poux_avant_tend, col =  "darkgreen")

Retranchement de la série (tendance)

Data_poux_sans_tend = zoom_Data_poux_ts - Data_poux_avant_tend

plot(Data_poux_sans_tend,
     xlab = "Années",
     ylab = "Indice")

Estimation de la saisonnalité à l’aide de moyennes saisonnières

# Calcul des moyennes saisonnières
moy_sais = aggregate(Data_poux_sans_tend~cycle(Data_poux_sans_tend), FUN=mean)

# Graphique 
plot(moy_sais,
     xlab = "Mois",
     ylab = "Indice",
     col = "darkblue",
     cex=2,
     pch=20,
     main = "Moyennes saisonnières")

Composante de saisonnalité à partir de ces moyennes

# Nombre de mois et d'années
S=12
nb_an = 9

# Définition des saisonnalités
poux_sais=rep(moy_sais$Data_poux_sans_tend, nb_an)
poux_sais = ts(poux_sais, start = c(2010,1), freq = 12)

# Graphique
plot(Data_poux_sans_tend,
     xlab="Années",
     ylab = " ",
     main = "Saisonnalité estimée")

lines(poux_sais, col = "darkgreen")

Retranchement de la saisonnalité à la série

poux_sans_tend_sais = Data_poux_sans_tend - poux_sais

plot(poux_sans_tend_sais,
     xlab="Années",
     ylab = "Indice",
     main = "Amelie Condette",
     col = "darkgreen")

Prévisions pour des valeurs futures de la série

# Prédicition de 4,5 ans
poux_pred=HoltWinters(zoom_Data_poux_ts)
pred_poux_holt=predict(poux_pred, n.ahead = 54, prediction.interval = TRUE)
plot(poux_pred, pred_poux_holt, main="Amelie Condette")

Erreurs de prévision

erreur_param = zoom_Data_poux_ts_apres - (pred_tend+pred_sais)
erreur_holt= zoom_Data_poux_ts_apres - pred_poux_holt
paste("Erreur de prévision moyenne du modèle paramétrique :", round(mean(erreur_param), 4), "Erreur de prévision moyenne du lissage de Holt-Winters :", round(mean(erreur_holt), 4))
## [1] "Erreur de prévision moyenne du modèle paramétrique : -0.0612 Erreur de prévision moyenne du lissage de Holt-Winters : 0.0161"