Avant de créer une chronique, il faut télécharger une librairie nommée zoo.

 

install.packages("zoo")

library(zoo)

 

 

Une chronique est un ensemble de couples valeur quantitative / date. Ce qui permet de voir une évolution des valeurs quantitatives en fonction du temps.

 

Chronique=zoo(Variable_Quantitative, order.by=Variable_Date) # On spécifie la variable de classe date avec order.by.

​

 

On peut visualiser une chronique à l’aide d’un graphique (plus précisément : un chronogramme). Pour cela, comme avec tous les graphiques, on utilise la fonction plot( ).

​

plot(Chronique, lwd=1, lty=1, col="blue", ylim=c(0,10000), xaxt="n") # ylim=c(0,10000) signifie que l’axe des ordonnées va de 0 à 10 000. xaxt="n" signifie que l’on va modifier l’axe des abscisses.

 

date = seq(as.Date("1982-01-01"),as.Date("2017-01-01"),by="5 years") # On crée une suite de dates.

axis(side=1,at=date,labels=format(date,"%b %Y")) # On applique la suite de dates à l’axe des abscisses.

 

 

On a donc l’évolution l’évolution des valeurs quantitatives entre janvier 1982 et janvier 2017.

Cependant, il arrive que l’on veuille proposer des prévisions de la chronique. Il est donc plus pertinent, de modéliser la chronique en sélectionnant uniquement les données les plus récentes.

On utilise donc la fonction window( ), pour spécifier une date de début. Puis on l’applique au graphique pour qu’il commence à cette date.

 

Chronique_Plus_Récente=window(Chronique,start=as.Date("2004-01-01")) # On commence la chronique au 1er janvier 2004.

​

plot(Chronique_Plus_Récente, lwd=1,lty=1,col="blue",ylim=c(2000,10000),xaxt="n")

date=seq(as.Date("2004-01-01"),as.Date("2017-01-01"),by="2 years")

axis(side=1,at=date,labels=format(date,"%b %Y")) # Même principe que précédemment. On modifie juste la date du début de la chronique.

 

 

Pour obtenir des caractéristiques numériques sur une série, on utilise les fonctions coredata( ) et index( ) permettant respectivement d’extraire les données brutes et les indices temporels de la chronique.

 

summary(coredata(Chronique_Plus_Récente)) # Caractéristiques sur les données brutes.

 

Moyenne_Par_Année=aggregate(Chronique_Plus_Récente, by=format(index(Chronique_Plus_Récente) ,"%Y"), FUN=mean) # Calcul de la moyenne des valeurs quantitatives par année.

 

 

On peut aussi représenter cette moyenne à l’aide d’un graphique.

 

plot(Moyenne_Par_Année) # Graphique des moyennes par année.

text(2008,5000,labels="Crise économique\n mondiale en 2008") # Texte dont les coordonnées sont 2008 en abscisse et 5 000 en ordonnée.

 

 

 

Pour obtenir une prévision de la chronique, il faut estimer la tendance, les saisonnalités et les résidus.

On va commencer par faire une estimation non paramétrique (lissage) de la tendance. Pour cela, on utilise la fonction loess( ).

 

time=as.numeric(index(Chronique_Plus_Récente)) # Indices temporels.

values=coredata(Chronique_Plus_Récente) # Données brutes.

smooth=loess(values~time,span=0.25,degree=1) # Estimation de la tendance en fonction des données brutes et des indices temporels de la chronique.

smooth=zoo(smooth$fitted,order.by=index(Chronique_Plus_Récente))

 

# Graphique.

plot(Chronique_Plus_Récente, col="blue", ylim=c(4000,8500), xaxt="n")

date = seq(as.Date("2004-01-01"),as.Date("2017-01-01"),by="2 years")

axis(side=1,at=date,labels=format(date,"%b %Y"))

 

# Estimation de la tendance par lissage.

lines(smooth,col="red",lwd=2)

​

# Création d’une légende.

legend("topleft",legend=c("Série brute","Série lissée"), col=c("blue","red"), lty=1, bty="n") 

​

​

On va maintenant modéliser la courbe de la tendance avec un modèle de régression linéaire par morceaux. On va donc définir des points de rupture (dans notre cas, on en a deux). Vous l'aurez compris une régression linéaire par morceaux est différent d'une courbe de lissage.

 

# Définition de deux points de rupture.

t1 = as.Date("2008-01-01") # Création de la date du premier point de rupture.

t2 = as.Date("2009-11-01") # Création de la date du deuxième point de rupture.

time = index(Chronique_Plus_Récente)

Ind1 = time > t1 ; time1 = (time-t1)*Ind1

Ind2 = time > t2 ; time2 = (time-t2)*Ind2

 

# Modélisation de la courbe de lisage par le modèle modèle de régression linéaire par morceaux.

values = coredata(smooth)

model1 = lm(values~time+time1+time2)

summary(model1)

trend = zoo(model1$fitted,order.by=index(Chronique_Plus_Récente))

 

# Graphique pour comparer l’estimation non paramétrique de la tendance (lissage) et le modèle de régression linéaire par morceaux.

plot(Chronique_Plus_Récente, lwd=1, lty=1, col="blue", ylim=c(4000,7000), xaxt="n", type="n")

​

date = seq(as.Date("2004-01-01"),as.Date("2017-01-01"),by="2 years")

axis(side=1,at=date,labels=format(date,"%b %Y"))

abline(v=date,col="gray50",lty=3)

lines(smooth,col="green4",lwd=2) # Estimation non paramétrique de la tendance.

lines(trend,col="red3",lwd=2) # Modèle de régression linéaire par morceaux.

legend("topleft",legend=c("Lissage","Modèle linéaire par morceaux"), col=c("green4","red3"), lty=1, bty="n") # Création d’une légende.

 

 

On a modélisé la composante tendance maintenant on va modélisé la composante saisonnière. Pour cela, vu que l’on est dans le cadre d’un schéma additif, on va retirer la composante tendancielle à la chronique.

 

# Graphique pour visualiser juste la composante saisonnière sans tenir compte la composante tendancielle.

Saison=Chronique_Plus_Récente - trend

plot(Saison,lwd=1,lty=1,col="blue",xaxt="n")

 

date = seq(as.Date("2004-01-01"),as.Date("2017-01-01"),by="2 years")

axis(side=1,at=date,labels=format(date,"%b %Y"))

abline(v=date,col="gray50",lty=3)

abline(h=0,lty=1)

 

​

Maintenant que l’on a vu qu’il y avait bien une saisonnalité, on va visualiser la distribution des données par mois et non plus par année.

 

# Visualisation des variations saisonnières.

Mois=aggregate(Saison, by=format(index(Chronique_Plus_Récente),"%m"), FUN=summary)

 

install.packages("lubridate")

library(lubridate) # Librairie qui permet d’utiliser la fonction month().

month = format(seq(as.Date("2016-01-01"),as.Date("2016-12-01"),by="month"),"%b")

boxplot(coredata(Saison)~month(Saison),col="orange",names=month)

points(1:12,Mois$Mean,type="b",col="green4",lwd=2)

abline(h=0,lty=3)

 

 

On va maintenant rechercher précisément la date à laquelle apparaissent les valeurs atypiques.

 

# Identification des valeurs atypiques.

which(Saison < -500 & month(Saison)==4)

Saison[76]

which( Saison < 0 & month(Saison)==9)

Saison[129]

 

Afin d’éviter que les valeurs atypiques n’ait un impact trop important sur l’estimation des coefficients saisonniers, on va remplacer ces valeurs par la médiane mensuelle correspondante.

 

# Traitement des valeurs atypiques.

Saison[76]=Mois[4,3]

Saison[129]=Mois[9,3]

 

 

On va maintenant déterminer une estimation des coefficients saisonniers à l’aide d’un modèle de régression linéaire.

 

# Estimation des coefficients saisonniers .

datamonth=rep(month,13)

datamonth=factor(datamonth,levels=month)

model2=lm(coredata(season2004)~-1+datamonth)

summary(model2)

 

 

On a obtenu des coefficients qualifiés de provisoires car ils ne vérifient par forcément l’hypothèse de conversation des aires. On va donc centrer les valeurs afin d’obtenir les coefficients saisonniers définitifs.

 

# Coefficients saisonniers définitifs.

Coef_Définitifs=model2$coef – mean(model2$coef)

 

​

On va maintenant visualiser l’ajustement saisonnier pour voir si les coefficients saisonniers définitifs sont bien ajustés.

 

season=zoo(rep(coefdef,13),order.by=index(Chronique_Plus_Récente))

plot(Saison,lwd=1,lty=1,col="blue",xaxt="n")

date=seq(as.Date("2004-01-01"), as.Date("2017-01-01"), by="2 years")

axis(side=1, at=date, labels=format(date,"%b %Y"))

abline(v=date, col="gray50", lty=3)

abline(h=0)

lines(season, col="red3", lwd=2, lty=4)

 

 

On va visualiser la série ajustée. Pour un schéma additif, cela correspond à la somme de la composante tendancielle et de la composante saisonière.

 

# Série ajustée.

adjusted = trend + season

plot(Chronique_Plus_Récente,lwd=1,lty=1,col="blue",ylim=c(3000,9000),xaxt="n")

date = seq(as.Date("2004-01-01"),as.Date("2017-01-01"),by="2 years")

axis(side=1,at=date,labels=format(date,"%b %Y"))

abline(v=date,col="gray50",lty=3)

lines(adjusted,col="red3",lwd=2)

legend("topleft",legend=c("Série brute","Série ajustée"), col=c("blue","red3"), lty=1, bty="n")

 

 

On a fait une estimation de la tendance puis de la saisonnalité, il reste plus qu’à faire une analyse de la série de résidus.

 

# Série des résidus et des résidus standardisés.

residual = Chronique_Plus_Récente - adjusted

residualstand = residual/sd(residual)

plot(residualstand,lwd=1.5,lty=1,col="blue",ylim=c(-6,+4),xaxt="n")

date = seq(as.Date("2004-01-01"),as.Date("2017-01-01"),by="2 years")

axis(side=1,at=date,labels=format(date,"%b %Y"))

abline(v=date,col="gray50",lty=3)

abline(h=0,lty=2)

abline(h=c(-2,+2),col="red3",lty=3)

 

 

Avant de poursuivre l’analyse des résidus, on va remplacer les données atypiques par l’indicateur de moyenne.

 

# Remplacement des données atypiques.

index = which(as.numeric(residualstand) < -3)

residual[index] = mean(residual[-index])

residualstand = residual/sd(residual)

 

​

On va réaliser un Quantile-Quantile Plot afin de confronter la série des résidus standardisés à un bruit blanc.

 

# Quantile-Quantile Plot.

install.packages("car")

library(car)

qqPlot(as.numeric(residualstand),pch=20,col="blue3", xlab="Quantiles théoriques", ylab="Quantiles empiriques")

 

 

Il peut être utile de réaliser un test de normalité de Shapiro-Wilk. Pour cela, on utilise la fonction shapiro.test( ).

 

data=as.numeric(residualstand)

shapiro.test(data)

 

 

 

On va maintenant représenter la fonction d’autocorrélation afin de la comparer à celle caractéristique d’un bruit blanc.

 

# Fonction d’autocorrélation.

acf(as.numeric(residualstand), lag.max=48, type="correlation", col="red4", lwd=2)

 

 

On va maintenant passer à la réaliser d’une prévision de la chronique sur les 4 prochaines années (pour janvier 2017 à décembre 2020).

 

# Préparations pour la prévision.

time=seq(as.Date("2017-01-01"), as.Date("2020-12-01"), by="months")

t1=as.Date("2008-01-01")

t2=as.Date("2009-11-01")

Ind1=time > t1 ; time1 = (time-t1)*Ind1

Ind2=time > t2 ; time2 = (time-t2)*Ind2

df=data.frame(time,time1,time2)

 

Il faut faire une extrapolation de la composante tendancielle et modifier l’extrapolation du trend. Pour cela, on fera des données du data.frame df et du modèle model1.

 

# Préparations pour la prévision.

predict.trend=predict(model1, newdata=df)

predict.trend=zoo(predict.trend, order.by=time)

predict=predict.trend + rep(coefdef,4)

print(round(predict,0))

 

Enfin, on va représenter graphiquement ces prévisions afin de les comparer à la chronique.

 

plot(Chronique_Plus_Récente, lwd=1, lty=1, col="blue", ylim=c(3000,10000), main="Fréquentation des passagers - Aéroports de Paris - Vols internationaux", xlab="Janvier 2004 - Décembre 2016", ylab="Nombre de passagers (en milliers)", xlim=c(as.Date("2004-01-01"),as.Date("2020-12-01")), xaxt="n")

 

date = seq(as.Date("2004-01-01"), as.Date("2020-12-01"), by="2 years")

axis(side=1, at=date, labels=format(date,"%b %Y"))

abline(v=date, col="gray50", lty=3)

lines(predict, col="red3", lwd=2, lty=3)

legend("topleft", legend=c("Chronique","Prévision"), col=c("blue","red3"), lty=c(1,3), bty="n") # Création de la légende.