# réglage du path: à adapter setwd("chemin/vers/votre/repertoire") # exercice 1: AirPassengers ########################### AirPassengers class(AirPassengers) start(AirPassengers) end(AirPassengers) frequency(AirPassengers) # frequency: nombre d'observation par unité de temps (12 mois dans l'année) plot(AirPassengers, ylab="Passengers (1000s)", type="o", pch =20) grid() # remarquer les maxima locaux en fin d'annee, au printemps (paques), et surtout en été # le modele multiplicatif X_t = T_t * S_t * u_t s'impose # car l'amplitude de la composante saisonniere ET des fluctuations aleatoires # augmentent avec la tendance # (cf graphique dans les slides) logAirPassengers <- log(AirPassengers) plot(logAirPassengers, ylab="Passengers (1000s) (log scale)", type="o", pch =20) grid() # après transformation, l'amplitude des comp. saisonnière et aléatoire # ne dépendent plus du niveau de y # le graphe ressemble davantage au graphe du modèle additif X_t = T_t + S_t + u_t # on verifie visuellement que le filtrage arithmetique donné dans le cours supprime la saisonnalité moymob <- filter (logAirPassengers, filter = c(1/24, rep(1/12,11), 1/24)) # graphes : plot(logAirPassengers, col="blue", ylab="Passengers (log scale)", type="o", pch=20) lines(moymob,col="red",type="o",pch=20) grid() # la suppression de la composante saisonnière permet de mettre en évidence un ralentissement du trafic en 1953 et 1957-1958 # ce qui n'était pas vraiment visible sur la chronique de départ windows() # créer un nouveau graphique; quartz() sous Mac et X11() sous Linux ? # calcul de la (pseudo) composante saisonnière saison <- logAirPassengers - moymob plot(saison, col="green", type="o",pch=20) grid() saison # à titre informatif (à faire dans l'exercice 3, où il y a les détails) # la transformée de Fourier étant linéaire, le périodogramme est à utiliser avec un modèle additif library(TSA) p <- periodogram(logAirPassengers) # on cherche la période correspondant au maximum dd <- data.frame(freq=p$freq, spec=p$spec, period=1/p$freq) dd <- dd[-c(1:4),] # on retranche les 4 premières valeurs order <- dd[order(-dd$spec),] top7 <- head(order, 7) top7 # pic "important" pour une période de 12 mois: c'est la période cherchée # décomposition decompAP <- decompose(AirPassengers,type="multiplicative") plot(decompAP) # la méthode plus sophistiquée (que notre moyenne mobile arithmétique) utilisée dans # decompose permet de trouver une composante saisonnière effectivement périodique # observer que comp. saisonnière et aléatoire sont centrées sur 1 (cas multiplicatif) # on retrouve les ralentissements dans la tendance en 1953 et 1958-1959. # reconstruction de chroniques qui nous intéressent: decompAP.sansbruit <- decompAP$trend*decompAP$seasonal decompAP.cvs <- decompAP$trend*decompAP$random plot(AirPassengers, col="blue", ylab="Passengers (log scale)", type="o", pch=20) lines(decompAP.sansbruit,col="red",type="o",pch=20) lines(decompAP.cvs,col="green",type="o",pch=20) # intérêt des données corrigées des variations saisonnières (CVS): # on voit par exemple des baisses ponctuelles (58-59, 60 par exemple) qui étaient masquées par la # composante saisonnière windows() plot(decompAP$seasonal) grid() # on voit mieux les maxima locaux en fin d'annee, au printemps (paques), et surtout en été # à comparer à: (par défaut, type="additive") decompAPadd <- decompose(AirPassengers) plot(decompAPadd) # random a un aspect assez "périodique", ce qui n'est pas conforme à nos hypothèses # remarquer en particulier que les composantes saisonnières et irrégulières # oscillent autour de 1 dans le modèle multiplicatif, # et autour de 0 dans le modèle additif # Exercice 2: IPI ################## IPI <- read.table(file="ipi.txt",sep='\t',header=TRUE) IPI <- ts(IPI, frequency = 4, start = c(1963, 1)) plot(IPI, ylab="Indice production industrielle", type="o", pch =20) grid() moymob <- filter (IPI, filter = c(1/8, rep(1/4,3), 1/8)) plot(IPI, col="blue", ylab="Indice production industrielle", type="o", pch=20) lines(moymob,col="red",type="o",pch=20) grid() # on voit déjà des ruptures de tendance sans doute causés par un phénomène exogène saison <- IPI - moymob plot(saison, col="green", type="o",pch=20) # la composante saisonnière n'est pas franchement périodique... # la représentation graphique initiale suggère plutôt un modèle multiplicatif decompIPI <- decompose(IPI,type="multiplicative") plot(decompIPI) # effet mai 1968 dans résidus + chocs pétroliers 1973 et 1979 avec effet retard dans tendance # (alors que ces effets n'étaient pas clairs dans la série "brute") #IPI.sansbruit <- decompIPI$trend*decompIPI$seasonal IPI.cvs <- decompIPI$trend*decompIPI$random windows() plot(IPI, col="blue", ylab="Indice production industrielle", type="o", pch=20) lines(IPI.cvs,col="red",type="o") grid() # on voit très bien dans les données CVS les effets mai 68 + chocs pétroliers # (bien mieux que dans la série brute) decompIPI$random # repérer 2ème trimestre 1968 qui est un point aberrant plot(decompIPI$random) grid() # Exercice 3: sunspot ##################### help(sunspot.month) # la chronique est mise à jour régulièrement sunspot.month # jusque 2013 sur la dernière version de R plot(sunspot.month) # il y a un visiblement une périodicité grid() library(TSA) # on voit la décroissance globale du "périodogramme" (ici plutôt le spectre de puissance) p <- periodogram(sunspot.month) # on voit un pic très prononcé signe de périodicité dd <- data.frame(freq=p$freq, spec=p$spec, period=1/p$freq) # on forme un dataframe dd <- dd[-c(1:4),] # on supprime les valeurs des 4 plus basses fréquences (ici ce n'est pas utile car le pic en indice 24 est largement plus grand) order <- dd[order(-dd$spec),] # on classe dans l'ordre décroissant des coefficients spectraux top10 <- head(order, 10) # on garde les 10 premières valeurs top10 # pic à 133.33 mois # durée du cycle en année: 133.33/12 # https://fr.wikipedia.org/wiki/Cycle_solaire ( 11.2 ans en moyenne)