1ier draft de la partie stat

module3/exo3/exercice_fr.Rmd

@@ -18,9 +18,11 @@ Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne co
 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)
@@ -48,7 +50,7 @@ plot(data_clean$DateExcel,data_clean$`CO2[ppm]`)
 
 
 
-Bon on va les détecter et les virer y'en a pas non plus trop donc autant faire les choses bien
+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,]
@@ -67,8 +69,76 @@ convert_iso_date = function(y,m) {
 
 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 je pense, sauf si vous en avez déjà marre mais ça serait dommage d'arrêter en si bon chemin. Cette étude est vraiment pleine de suspens...
 
 
-## Modèle statistique
\ No newline at end of file
+## Modèle statistique
+
+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 lisasge par noyaux mais bon.
+
+```{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')
+
+
+```
+
+
+
+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)+I(Yr^3),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')
+```
+
+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')
+
+```
+
+Sinon on aurait pu se concentrer sur la distribution des évolutions annuelles
+```{r}
+evol_ann<-diff(MA$MA)/MA$MA[-length(MA)]
+plot(MA$Yr[-1],evol_ann)
+plot(MA$Yr[-1],diff(MA$MA))
+
+
+reg_evol<-lm(evol_ann~MA$Yr[-1])
+summary(reg_evol)
+plot(reg_evol)
+
+
+```
\ No newline at end of file