Dans ce document nous allons présenter et analyser des données d'évolution de la concentration en CO2 atmosphérique provenant
du site de l'observatoire de Mauna Loa situé sur l'île d'Hawaï.
Dans ce document nous allons présenter et analyser des données d'évolution de la concentration en CO2 atmosphérique provenant du site de l'observatoire de Mauna Loa situé sur l'île d'Hawaï.
*C. D. Keeling, S. C. Piper, R. B. Bacastow, M. Wahlen, T. P. Whorf, M. Heimann, and
H. A. Meijer, Exchanges of atmospheric CO2 and 13CO2 with the terrestrial biosphere and
oceans from 1978 to 2000. I. Global aspects, SIO Reference Series, No. 01-06, Scripps
Institution of Oceanography, San Diego, 88 pages, 2001.*
*C. D. Keeling, S. C. Piper, R. B. Bacastow, M. Wahlen, T. P. Whorf, M. Heimann, and H. A. Meijer, Exchanges of atmospheric CO2 and 13CO2 with the terrestrial biosphere and oceans from 1978 to 2000. I. Global aspects, SIO Reference Series, No. 01-06, Scripps Institution of Oceanography, San Diego, 88 pages, 2001.*
## Chargement des données
...
...
@@ -63,37 +62,37 @@ data[na_records,]
**Premier problème rencontré** : les colonnes 3 et 4 sont des dates mais le format de la date n'est pas indiqué.
Après recherche voici la description des quatre première colonne :
- colonne 1 : année ;
- colonne 2 : mois ;
- colonne 3 : date en nombre de jours depuis une origine (comme dans Excel). Date d'origine 1899-12-30 ;
- colonne 4 : date en année décimal.
- colonne 1 : année ;
- colonne 2 : mois ;
- colonne 3 : date en nombre de jours depuis une origine (comme dans Excel). Date d'origine 1899-12-30 ;
- colonne 4 : date en année décimal.
## Conversion des colonnes en date
### Conversion des colonnes en date
Nous allons convertir la colonne 3 en date en indiquant la date 1899-12-30 comme origine.
Nous allons convertir la colonne 3 en date en indiquant la date 1899-12-30 comme origine.
Cette conversion renvoie le 15 du mois ce qui est conforme à la description dans le fichier "The monthly values have been adjusted to 24:00 hours on the 15th of each month"
Cette conversion renvoie le 15 du mois ce qui est conforme à la description dans le fichier "The monthly values have been adjusted to 24:00 hours on the 15th of each month".
Nous n'avons pas réussi à faire de même avec la colonne 4.
## Représentations graphique
## Représentations graphiques des données
### Réprésentation de l'évolution de la concentration de CO2 atmosphérique en fonction du temps
D'après la description des données contenu dans le fichier csv la colonne 5 CO2 correspond à la concentration
en CO2 : "Column 5 below gives monthly Mauna Loa CO2 concentrations in micro-mol CO2 per mole (ppm)".
D'après la description des données contenu dans le fichier csv la colonne 5 CO2 correspond à la concentration en CO2 : "Column 5 below gives monthly Mauna Loa CO2 concentrations in micro-mol CO2 per mole (ppm)".
Ce ne sont pas des données brutes mais des données recalculées à une échelle mensuelle.
mtext("Mauna Loa Observatory, Hawaii", side=3,line=2,cex=1.25,adj=0)
mtext("Monthly average carbone dioxide concentration", side=3,line=1,cex=1.25,adj=0)
mtext("Data from Scripps CO2 Program", side=3,lin=0,cex=.75,adj=0)
...
...
@@ -107,13 +106,15 @@ Nous voyons clairement deux structures dans les données, une forte tendance à
### Transformation des données en série temporelle
Pour caractériser ces deux composantes, nous allons considérer les données comme une série temporelle.
Nous allons transformer les données en un objet ts. Pour cela nous devons utiliser une série de données sans valeurs manquantes.
Plutôt que la colonne 5 nous allons donc utiliser la colonne 9 "CO2.1" qui correspond aux mêmes données que la colonne 5
mais où les données manquantes ont été remplacées en utilisant une "smoothed version of the data generated from a stiff cubic spline function plus 4-harmonic functions with linear gain". Cette procédure est décrite en détail dans le rapport cité en introduction de ce document.
Nous pouvons le vérifier en superposant les données remplies et les données avec données manquantes :
Nous allons transformer les données en un objet ts.
Pour cela nous devons utiliser une série de données sans valeurs manquantes.
Plutôt que la colonne 5 nous allons donc utiliser la colonne 9 "CO2.1" qui correspond aux mêmes données que la colonne 5 mais où les données manquantes ont été remplacées en utilisant une "smoothed version of the data generated from a stiff cubic spline function plus 4-harmonic functions with linear gain".
Cette procédure est décrite en détail dans le rapport cité en introduction de ce document.
Nous pouvons le vérifier en superposant les données remplies et les données avec données manquantes :
mtext("Mauna Loa Observatory, Hawaii", side=3,line=2,cex=1.25,adj=0)
mtext("Monthly average carbone dioxide concentration", side=3,line=1,cex=1.25,adj=0)
mtext("Data from Scripps CO2 Program", side=3,lin=0,cex=.75,adj=0)
...
...
@@ -124,7 +125,8 @@ legend("topleft",c("CO2 data with NA","CO2 data filled"),lty=1,col=c(1,2),bty="n
Nous constatons qu'il y a très peu de données manquantes et qu'elles se situent toutes entre 1958 et 1964 (nous les avions déjà identifiées plus haut).
legend("topleft",c("CO2 data with NA","CO2 data filled"),lty=1,col=c(1,2),bty="n")
```
Il nous faut également supprimer les premières et dernières données (manquantes)
allant de janvier 1958 à février 1958 puis de juin 2023 à décembre 2023.
Il nous faut également supprimer les premières et dernières données (manquantes) allant de janvier 1958 à février 1958 puis de juin 2023 à décembre 2023.
```{r}
data.trunc<-data[!is.na(data$CO2.1),]
head(data.trunc[,c("Date","CO2","CO2.1")],10)
tail(data.trunc[,c("Date","CO2","CO2.1")],10)
```
Nous pouvons maintenant transformer la colonne 9 en un objet ts.
Nous commençons par trier les données pour être sûr que les lignes sont bien ordonnées selon la date.
La série temporelle correspond à des données mensuelles sur plusieurs années.
La première donnée est le mois de mars 1958. La fin est le mois de mai 2023.
La première donnée est le mois de mars 1958.
La fin est le mois de mai 2023.
Nous créons notre objet ts en indiquant le début et la fréquence de la série.
```{r}
...
...
@@ -154,8 +157,10 @@ plot(data.ts)
```
### Décomposition de la série temporelle
La fonction stl du package stats permet de décomposer la série temporelle en utilisant un modèle loess pour ajuster la tendance.
Nous avons fixé le paramètre s.window à 13 ce qui signifie que la fonction utilise une fenêtre de 13 annés consécutives pour estimer chaque valeur de la composante saisonnière. Nous avons préféré cela plutôt que la valeur par défaut ("periodic") qui considère que la composante saisonnière ne varie pas au cours du temps.
La fonction stl() du package stats permet de décomposer la série temporelle en utilisant un modèle loess pour ajuster la tendance.
Nous avons fixé le paramètre s.window à 13 ce qui signifie que la fonction utilise une fenêtre de 13 annés consécutives pour estimer chaque valeur de la composante saisonnière.
Nous avons préféré cela plutôt que la valeur par défaut ("periodic") qui considère que la composante saisonnière ne varie pas au cours du temps.
La fonction renvoie la composante saisonnière, la tendance et les résidus.
```{r}
...
...
@@ -182,7 +187,8 @@ mtext("hors tendance",line=0)
```
Nous pouvons également donner quelques statistiques de cette composante saisonnières.
Dans une année le CO2 atmosphérique peut varier de `r round(min(data.trunc$seasonSTL),2)` ppm à `r round(max(data.trunc$seasonSTL),2)` ppm par rapport à la moyenne annuelle. L'écart-type de la composante saisonnière est de `r round(sd(data.trunc$seasonSTL),2)` ppm.
Dans une année le CO2 atmosphérique peut varier de `r round(min(data.trunc$seasonSTL),2)` ppm à `r round(max(data.trunc$seasonSTL),2)` ppm par rapport à la moyenne annuelle.
L'écart-type de la composante saisonnière est de `r round(sd(data.trunc$seasonSTL),2)` ppm.
Nous voyons sur le graphe précédent que cette variation intra-annuelle augmente au cours du temps.
La variance de la composante saisonnière n'est pas constante au cours du temps.
...
...
@@ -191,13 +197,15 @@ La variance de la composante saisonnière n'est pas constante au cours du temps.
Nous pouvons superposer la tendance aux données brutes ce qui correspond au même type de données que la colonne 6 seasonally : "Column 6 gives the same data after a seasonal adjustment to remove the quasi-regular seasonal cycle"
ce modèle est simpliste car par exemple nous n'avons pas tenu compte des observations répétées par années, c'est à dire de la présence de pseudoréplications, ni de l'autocorrélation temporelle.
Néanmoins, même s'il n'est pas parfait, ce modèle très simple semble bien convenir aux données.
### Prédiction jusqu'en 2025 à l'aide du modèle simple
Pour prédire la concentration de CO2 en 2025, nous allons utiliser le modèle précédent et l'appliquer jusqu'en 2025.
Au delà de mai 2023, nous n'avons plus de données.
Il s'agira d'extrapolations à l'aide du modèle simple.
Nous utilisons la fonction predict.lm qui permet d'obtenir l'intervalle de confiance des prédictions.
Nous commençons par construire un data.frame avec une colonne "time" sur l'ensemble de la période de mars 1958 à décembre 2025.
Puis nous utilisons ce data.frame pour prédire les valeurs de CO2 à l'aide du modèle simple.
Nous créons un data.frame spécifiquement pour les données au-delà de mai 2023 qui correspond aux valeurs prédites extrapolées.
Le graphe suivant représente la tendance "observée" sur la période connue, les valeurs de la tendance obtenues par ajustement du modèle simple et enfin les valeurs de la tendance prédite au-delà de mai 2023.
legend("topleft",c("tendance observée","valeurs prédites","intervalle de confiance des prédictions"),col=c("red","black","black"),lty=c(1,1,2),bty="n",cex=0.75)
**D'après le modèle simple, en première approximation la valeur de concentration moyenne en CO2 en 2025 serait de `r round(mean(myNewData2025$fit),2)` ppm plus ou moins `r round(mean(myNewData2025$fit)-mean(myNewData2025$lwr),2)` ppm.**
### Modélisation à l'aide d'un modèle stochastique
Pour aller un peu plus loin nous allons tenter d'utiliser un modèle stochastique de type SARIMA pour faire les prédictions.
En effet, nous avons vu précédemment que les données présentaient à la fois une tendance et une saisonnalité.
Nous pouvons également vérifier que les résidus présentent une autocorrélation et une autocorrélation partielle à l'aide des fonctions acf et pcaf.
Pour simplifier le processus de construction du modèle nous allons utiliser la fonction auto.arima() du package forecast qui permet de choisir le type de modèle approprié et d'estimer les différents paramètres du modèle stochastique.
```{r}
#install.packages("forecast")
library(forecast)
modelSarima<-auto.arima(data.ts,stepwise=FALSE)
summary(modelSarima)
```
Nous pouvons vérifier qu'il n'y a plus d'autocorrélation dans les résidus du modèle.
```{r}
par(mfrow=c(1,2))
acf(resid(modelSarima))
pacf(resid(modelSarima))
```
Pour prédire la concentration de CO2 jusqu'en 2025 nous allons utiliser la fonction forecast() du package forecast.
Pour aller jusqu'à décembre 2025 nous allons prédire sur 31 valeur supplémentaires (de juin 2023 à décembre 2025).
Nous pouvons représenter les prédictions à l'aide de la fonction plot.forecast() du package forecast qui représente les prédictions avec les intervalles de confiance à 80 et 95%.
```{r}
plot(forecast(modelSarima,h=31),include=50)
```
De la même façon que pour le modèle simple nous pouvons estimer la concentration de CO2 en 2025.
Nous récupérons les dates dans le data frame des prédictions du modèle SARIMA et nous créons un data.frame uniquement pour 2025.
**D'après le modèle SARIMA, la valeur de concentration moyenne en CO2 en 2025 serait de `r round(mean(predictionsSARIMA2025$Point.Forecast),2)` ppm plus ou moins `r round(mean(predictionsSARIMA2025$Point.Forecast)-mean(predictionsSARIMA2025$Lo.95),2)` ppm.**
Le modèle SARIMA prédit donc une concentration de CO2 très proche par rapport au modèle simple, très légèrement inférieure.
L'intervalle de confiance des prédictions est légèrement plus large avec le modèle SARIMA.
