---
title: 'Sujet 1 : Concentration de CO2 dans l''atmosphère depuis 1958'
author: "E. Guiffart"
date: "1 Juin 2020"
output:
  pdf_document: default
  pdf: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Préparation des données

Les données de *la concentration du CO2 dans l'atmosphère* sont disponibles sur le site web du [Scripps CO2 Program](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html), un site très intéressant comme on devrait en voir plus souvent de nos jours, sauf qu'il est américain. 

Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne correspond à une semaine de la période demandée. Nous téléchargeons toujours le jeu de données complet, qui commence en 1991 et se termine avec une semaine récente.

En pratique, nous chargons les données grâce à une copie du fichier enregistrée localement. Un code télécharge directement et enregistre une copie au cas où ce fichier n'exste pas encore en local.
```{r}
library(data.table)
library(timeSeries)
library(parsedate)


if (!file.exists("monthly_in_situ_co2_mlo.csv")) download.file("https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/monthly/monthly_in_situ_co2_mlo.csv", "monthly_in_situ_co2_mlo.csv", method="auto") 

data_url <- fread(file = "monthly_in_situ_co2_mlo.csv",skip=54, sep=",",strip.white=TRUE, header = FALSE)
```

Allez on visualise les données:
```{r}
head(data_url)
tail(data_url)
```

## Nettoyage des données
Va falloir nettoyer tout ça, notamment le nom des colonnes:
```{r}
new_names<-paste0(data_url[1,],data_url[2,],data_url[3,])
setnames(data_url,names(data_url),new_names)
data_clean<-data_url[-c(1,2,3),]
head(data_clean)
```

Bon ok c'est un peu mieux, moins confusant c'est déjà ça. Les 2 premières lignes sont suspectes. Voyons s'il ya en a d'autres. Pour ça, faisons un petit graph :
```{r}
plot(data_clean$DateExcel,data_clean$`CO2[ppm]`)
```



Bon on va les détecter et les modifier en les supprimant. 
```{r}
na_records = apply(data_clean, 1, function (x) any(is.na(x)+(x<0)))
data_clean[na_records,]

data_clean2<-data_clean[!na_records,]
plot(data_clean2$DateExcel,data_clean2$`CO2[ppm]`)

```

ok on est sur la bonne voie, on va s'en sortir. Reste les dates à mettre en meilleur format:
```{r}
convert_iso_date = function(y,m) {
	iso = paste0(y, "-", m, "-","01")
	as.character(parse_iso_8601(iso))
}

data_clean2$dateStd<-as.Date(convert_iso_date(data_clean2$Yr,data_clean2$Mn))
head(data_clean2)

plot(data_clean2$dateStd,data_clean2$`CO2[ppm]`)
```
On va pouvoir passer aux stats maintenant. Cette étude est vraiment pleine de suspens...


## Visualisation des oscillations périodiques et tendance

Regardons de plus près l'évolution de la série en focalisant sur des années, par exemple de 1980 à 1985:
```{r}


plot(data_clean2[Yr>1980 & Yr<1986,]$dateStd,data_clean2[Yr>1980 & Yr<1986,]$`CO2[ppm]`)
```

Nous observons bien une saisonalité de période annuelle. On ne va que s'intéresser dans la suite de cette étude aux données brutes c'est à dire le champ "*CO2(pm)*" car je ne sais pas trop à quoi correpond le reste des champs. Pour isoler la partie tendance, nous allons utiliser une technique simple de moyenne mobile. C'est un estimateur non-paramétrique de la tendance, au sens ou nous ne supposons pas de structure a-priori de cette tendance (par ex. linéaire ou polynomiale). On aurait aussi pu utiliser un lissage par noyaux mais bon restons sur des choses basiques.

```{r}

MA<-data_clean2[, .(MA=mean(as.numeric(`CO2[ppm]`))),by=Yr]
MA[,Yr_date:=as.Date(convert_iso_date(Yr,"06"))]
MA$Yr<-as.numeric(MA$Yr)
tail(MA)

plot(data_clean2$dateStd,data_clean2$`CO2[ppm]`,type='l')
lines(MA$Yr_date,MA$MA,col='red')


```


Nous obtenons la tendance de la série corrigée des effets de saisonalité

## Modèle statistique et calibrage de la tendance (contribution lente)
Nous allons construire un modèle simple de la tendance qui servira à prédire le niveau de CO2 dans l'atmosphère pour les années futures.
Faisons une régression sur la tendance. On pourrait enlever les premiers et derniers points car ce ne sont pas des donnés complètes.
```{r}

plot(MA$Yr,log(MA$MA),col='red',type='l')
reg_lin_log<-lm(MA~log(Yr),data=MA)
summary(reg_lin_log)
plot(reg_lin_log)

reg_lin<-lm(MA~Yr+I(Yr^2),data=MA)
summary(reg_lin)
plot(reg_lin)
tail(predict(reg_lin))
plot(MA$Yr,MA$MA,type='l')
lines(MA$Yr,predict(reg_lin),col='red')
```

Le modèle de régression est donc créé.


## Prédiction / Extrapolation
Il est temps de passser à l'extrapolation:
```{r}

dates_futur<-as.Date(convert_iso_date(seq(2021,2025),"06"))
futur<-data.table(Yr=seq(2021,2025))

pred<-predict(reg_lin,futur,interval = "prediction")
plot(MA$Yr,MA$MA,type='l',xlim=c(1960,2025),ylim=c(310,450))
lines(seq(2021,2025),pred[,1],col='red')

```
La courbe en rouge représente les valeurs futures jusqu'à 2025.