Sujet

En 1958, Charles David Keeling a initié une mesure de la concentration de CO2 dans l’atmosphère à l’observatoire de Mauna Loa, Hawaii, États-Unis qui continue jusqu’à aujourd’hui. L’objectif initial était d’étudier la variation saisonnière, mais l’intérêt s’est déplacé plus tard vers l’étude de la tendance croissante dans le contexte du changement climatique. En honneur à Keeling, ce jeu de données est souvent appelé “Keeling Curve” (voir (https://en.wikipedia.org/wiki/Keeling_Curve) pour l’histoire et l’importance de ces données).

Les données sont disponibles sur le site Web de l’ Institut Scripps . J’ai utilisé par erreur les observations mensuelles (je n’ai pas le temps de refaire tout le travail en utilisant les données hebdomadaires, ce sont d’ailleurs les données mensuelles qui sont utilisées en général dans les études d’impact climatiques ) observations mensuelles. Attention, ce fichier est mis à jour régulièrement avec de nouvelles observations. Notez donc bien la date du téléchargement, et gardez une copie locale de la version précise que vous analysez. Faites aussi attention aux données manquantes.

Pré-requis

Traitement de suites chronologiques

Quelques références:

(https://github.com/Peymankor/Data-Science_Portfolio/blob/master/Time%20Series%20Analysis-Historical%20Co2/medium_post.Rmd)

Pré-traitements sur les données

Extrait des métadonnées dans le header du fichier csv téléchargé.

The data file below contains 10 columns.

Missing values are denoted by -99.99

CO2 concentrations are measured on the ‘08A’ calibration scale

Chargement des librairies R, du fichier de données, renommage des colonnes, premières statistiques

setwd("C:/Users/hraynal/EspaceTravailBadet/FORMATION_PERSO/2020-Reproductibilite/mooc-rr/module3/exo3")
library(tidyverse)
library(forecast)
library(lubridate)
library(car)
library(scales)
library(patchwork)
library(kableExtra)

dataCO2 <- read.csv("monthly_in_situ_co2_mlo.csv", sep="," ,skip = 57)
colnames(dataCO2) <- c("Year", "Month","Date1", "Date2", "ObsCO2", "SeasAdjCO2","SplineAdjCO2", "SplineAdjCO2Trend", "ObsCO2Comp", "SeasAdjCO2Comp")
summary(dataCO2)
##       Year          Month            Date1           Date2     
##  Min.   :1958   Min.   : 1.000   Min.   :21231   Min.   :1958  
##  1st Qu.:1973   1st Qu.: 4.000   1st Qu.:26968   1st Qu.:1974  
##  Median :1989   Median : 7.000   Median :32704   Median :1990  
##  Mean   :1989   Mean   : 6.507   Mean   :32705   Mean   :1990  
##  3rd Qu.:2005   3rd Qu.: 9.500   3rd Qu.:38442   3rd Qu.:2005  
##  Max.   :2020   Max.   :12.000   Max.   :44180   Max.   :2021  
##      ObsCO2         SeasAdjCO2      SplineAdjCO2    SplineAdjCO2Trend
##  Min.   :-99.99   Min.   :-99.99   Min.   :-99.99   Min.   :-99.99   
##  1st Qu.:328.40   1st Qu.:328.70   1st Qu.:328.46   1st Qu.:328.82   
##  Median :351.34   Median :352.13   Median :351.33   Median :352.03   
##  Mean   :346.18   Mean   :346.18   Mean   :348.95   Mean   :348.95   
##  3rd Qu.:377.55   3rd Qu.:377.35   3rd Qu.:377.69   3rd Qu.:377.37   
##  Max.   :414.83   Max.   :413.33   Max.   :414.94   Max.   :413.35   
##    ObsCO2Comp     SeasAdjCO2Comp  
##  Min.   :-99.99   Min.   :-99.99  
##  1st Qu.:328.40   1st Qu.:328.70  
##  Median :351.34   Median :352.13  
##  Mean   :348.96   Mean   :348.95  
##  3rd Qu.:377.55   3rd Qu.:377.35  
##  Max.   :414.83   Max.   :413.33
dataCO2$Date <- ymd(paste0(dataCO2$Year, " ", dataCO2$Month, " ", "15"))

Traitement des valeurs manquantes notées -99.99

Il s’avère que dans la série des valeurs observées complétée par les valeurs interpolées, il reste des observations avec des valeurs manquantes (-99.99). On supprime ces observations.

dataCO2 <- dataCO2[dataCO2$ObsCO2Comp != "-99.99", ]

Création d’une variable date au format YYYY MM DD

dataCO2$Date <- ymd(paste0(dataCO2$Year, "-", dataCO2$Month, "-", "15"))

Représentations graphiques des résultats

ggplot(dataCO2,aes(Date, dataCO2$ObsCO2Comp)) +
    geom_line(color='orange') +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "5 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1),axis.title.y = element_text(size = 10)) +
  ggtitle("Graphique 1: Evolution au cours du temps de la teneur en CO2")

dataCO2_by_year <- dataCO2 %>% group_by("Year")
ggplot(dataCO2_by_year, aes(dataCO2_by_year$Month,dataCO2_by_year$ObsCO2Comp )) +
         geom_line(aes( group = dataCO2_by_year$Year , colour=dataCO2_by_year$Year)) +
         xlab("Month")+
         ylab("CO2 Concentration (ppm)") +
         ggtitle("Graphique 2: Evolution saisonnière de la concentration en CO2, par année")

- On déduit du graphique 1 que la série n’est pas stationnaire, elle présente une tendance (augmentation du taux de CO2 au cours du temps).

Question 1: Réalisez un graphique qui vous montrera une oscillation périodique superposée à une évolution systématique plus lente.

# Définition en tant que série temporelle
Co2_train <- ts(dataCO2$ObsCO2Comp, start = c(1958,3), frequency = 12)
T <- time(Co2_train)
base <- data.frame(Y=dataCO2$ObsCO2Comp,X1=as.numeric(T))

#essai d'ajustement à un modèle de régression linéaire
reg1 <- lm(Y~X1,data=base)
tendance1 <-predict(reg1)
T <-as.numeric(T)

#Calcul des moyennes annuelles de teneur en CO2
# <- cbind(  levels(factor(dataCO2$Year)),  tapply(dataCO2$ObsCO2Comp, dataCO2$Year, mean)) 

MeanCO2 <- tapply(dataCO2$ObsCO2Comp, dataCO2$Year, mean)

tmp <- ymd(paste0(unique(dataCO2$Year), "-06-15")) 
MeanCO2Year <- data.frame(tmp , as.numeric(MeanCO2) )
colnames(MeanCO2Year) <- c("Year", "MeanCO2") 
dataCO2M <- merge( dataCO2, MeanCO2Year, by="Year")



ggplot() +
   geom_line(aes( x=dataCO2$Date, y=dataCO2$ObsCO2Comp), color='orange') +
   geom_line(aes(x=MeanCO2Year$Year, y=MeanCO2Year$MeanCO2  ), color="blue") +
  geom_line(aes(x=dataCO2$Date, y=tendance1  ), color="grey") +
   #geom_line(aes(x=T, y=tendance1  ), color="grey") +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "5 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1),axis.title.y = element_text(size = 10)) +
  ggtitle("Graphique 1: Evolution au cours du temps de la teneur en CO2")

Modélisation de la série temporelle

On va chercher à modéliser la série en utilisant une catégorie de modèles qui cherche à déterminer chaque valeur de la série en fonction des valeurs qui la précède (yt = f(yt-1, yt-2, …)). C’est le cas des modèles ARIMA (“AutoRegressive – Integrated – Moving Average”). Cette catégorie de modèles a été popularisée et formalisée par Box et Jenkins (1976).

Les processus autorégressifs supposent que chaque point peut être prédit par la somme pondérée d’un ensemble de points précédents, plus un terme aléatoire d’erreur.

Le processus d’intégration suppose que chaque point présente une différence constante avec le point précédent.

Les processus de moyenne mobile supposent que chaque point est fonction des erreurs entachant les points précédant, plus sa propre erreur.

Un modèle ARIMA est étiqueté comme modèle ARIMA (p,d,q), dans lequel:

On va utiliser la fonction auto.arima() qui est un algorithme automatique qui permet de trouver le modèle ARIMA qui estime le mieux la série temporelle.

Co2_train <- ts(dataCO2$ObsCO2Comp, start = c(1958,3), frequency = 12)
Co2_train %>% ggtsdisplay()

auto.arima(Co2_train)
## Series: Co2_train 
## ARIMA(1,1,1)(2,1,2)[12] 
## 
## Coefficients:
##          ar1      ma1     sar1     sar2    sma1     sma2
##       0.2188  -0.5797  -0.3038  -0.0090  -0.545  -0.2643
## s.e.  0.0891   0.0746      NaN   0.0408     NaN      NaN
## 
## sigma^2 estimated as 0.09695:  log likelihood=-182.81
## AIC=379.62   AICc=379.77   BIC=411.79

En se basant sur le résultat de l’analyse précédente, on retient le modèle \(ARIMA(1,1,1)(2,1,2)[12]\). On peut ensuite vérifier que les résidus ressemblent à un bruit blanc ou non.

(fit_minaicc <- Arima(Co2_train, order=c(1,1,1),seasonal=list(order=c(2,1,2),period=12),
                  lambda = "auto"
              ))
## Series: Co2_train 
## ARIMA(1,1,1)(2,1,2)[12] 
## Box Cox transformation: lambda= 0.3024523 
## 
## Coefficients:
##          ar1      ma1     sar1     sar2     sma1     sma2
##       0.2105  -0.5733  -0.3345  -0.0088  -0.5363  -0.3081
## s.e.  0.0897   0.0757      NaN      NaN      NaN      NaN
## 
## sigma^2 estimated as 2.992e-05:  log likelihood=2824.66
## AIC=-5635.33   AICc=-5635.17   BIC=-5603.16
checkresiduals(fit_minaicc, lag=36)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(2,1,2)[12]
## Q* = 25.374, df = 30, p-value = 0.7067
## 
## Model df: 6.   Total lags used: 36
fit_minaicc$aicc
## [1] -5635.174

Les graphiques précédents montre que le résidu ressemble suffisamment au bruit blanc, que la valeur p est élevée et que le modèle a passé le test de la Ljong-Box.

autoplot(forecast(fit_minaicc, h=200))