---
title : Analyse de la concentration de CO2 dans l'atmosphère depuis 1958
author : François Févotte
date: avril 2020
options:
  css: skeleton_css.css
  template: julia_html.tpl
---

Ce document a été réalisé à l'occasion du [MOOC "Recherche
Reproductible"](https://www.fun-mooc.fr/courses/course-v1:inria+41016+self-paced/info). Il
a donc plus vocation à fournir des résultats *reproductibles* que des résultats
*intéressants*. N'étant pas un scientifique des données, il est tout à fait
possible que l'analyse réalisée ici soit caduque ; je plaide pour l'indulgence
du lecteur !

Ce document a été automatiquement généré à partir du fichier
[`exercice.jmd`](exercice.jmd), que le lecteur intéressé pourra consulter
directement. Une [annexe](repro.html) détaille toutes les étapes nécessaires
afin de reproduire l'analyse et re-générer intégralement le présent
document. Cette annexe, nommée [`repro.html`](repro.html), devrait se trouver
dans le même dépôt git et même répertoire que le présent document.

------

**Résumé:**

Depuis 1958, la concentration atmosphérique en CO2 est mesurée régulièrement à
l'observatoire de Mauna Loa, à Hawaï. Ces mesures, initiées par [Charles
D. Keeling](https://en.wikipedia.org/wiki/Charles_David_Keeling), ont été les
parmi les premières à mettre en évidence l'accroissement rapide de la
concentration en CO2 dans l'atmosphère, un constat qui a par la suite permis
d'attirer l'attention de la communauté scientifique (puis par la suite du grand
public) sur l'impact des activités humaines sur l'atmosphère.

![Courbe de
Keeling](https://upload.wikimedia.org/wikipedia/commons/thumb/c/c5/Mauna_Loa_CO2_monthly_mean_concentration.svg/480px-Mauna_Loa_CO2_monthly_mean_concentration.svg.png)

```julia; echo=false;
monroe2015() = print("[[monroe2015](https://scripps.ucsd.edu/programs/keelingcurve/2015/02/12/is-the-rate-of-co2-growth-slowing-or-speeding-up/)]");
```

La [courbe de Keeling](https://en.wikipedia.org/wiki/Keeling_Curve) fait
apparaître une augmentation de la concentration atmosphérique de l'ordre de 1 à
2 ppm/an `j monroe2015()`. Elle fait aussi apparaître une variation
cyclique d'environ 5 ppm chaque année, qui rend l'analyse de la tendance
systématique haussière plus délicate.

Nous nous proposons ici de refaire pas à pas l'étude de ce jeu de données. Nous
séparons la composante cyclique de la composante tendancielle, afin de
caractériser cette dernière. Les paramètres de notre modèle correspondent à un
taux d'augmentation du CO2 atmosphérique de 0.8 ppm/an en 1959, et 1.5 ppm/an en
2015 (cette dernière valeur étant un peu sous-estimée par rapport aux 2.25
ppm/an annoncés dans `j monroe2015()`).

Notre modèle est construit en utilisant les données de 1959 à 2014. Il est
validé par comparaison aux données de 2015 à nos jours. Une extrapolation sur la
période 2020-2025 est aussi réalisée afin de permettre une comparaison à de
futures mesures.


# Gestion des dépendances


## Environnement

```julia; echo=false; results="hidden"
@info "* Instantiating project"
```

Nous utilisons Julia dans sa version 1.4.0, sur une architecture matérielle de
type x86 (64 bits)
```julia
using InteractiveUtils
versioninfo()
```

Il nous faut maintenant charger l'environnement logiciel de notre étude, c'est à
dire toutes les bibliothèques sur lesquelles elle s'appuie, dans les bonnes
versions. La liste des dépendances directes est contenue dans le fichier
`Project.toml`, et complétée par le fichier `Manifest.toml` (qui liste toutes
les dépendances directes et indirectes, accompagnées de leurs numéros de version
précis).

Toutes ces informations permettent au gestionnaire de paquets de re-créer un
environnement logiciel identique à celui qui a été utilisé pour développer cette
analyse.
```julia
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.status()
```


## Chargement des dépendances

```julia; echo=false; results="hidden"
@info "* Loading dependencies (this might take a while...)"
```

Tant que nous y sommes, profitons en pour charger dès maintenant les paquets
dont nous aurons besoin.
```julia
import HTTP
import CSV
import DataInterpolations; const DI=DataInterpolations
using GLM
using Printf
using Dates
using DataFrames
using Statistics
using Plots; plotly()
```

```julia; echo=false; results="hidden"
include("utils.jl")
```

# Données d'entrée

Nos données d'entrée proviennent du programme [Scripps
CO2](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html). Nous
fondons l'analyse sur le jeu de données contenant des observations hebdomadaires.


## Téléchargement

```julia; echo=false; results="hidden"
@info "* Retrieving data"
```

Le jeu de données est téléchargé une seule fois ; c'est une copie locale qui
sert à réaliser l'analyse. Ceci permet de garantir la reproductibilité des
données utilisées pour l'analyse, qui seront stockées dans git aux côtés du
présent document.

Il est toutefois possible de forcer le téléchargement en positionnant la
variable `force_download=true`. Ceci permet notamment d'actualiser le jeu de
données.

```julia; results="hidden"
const data_url = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv"
const data_file = "weekly_in_situ_co2_mlo.csv"
const force_download = false
```

```julia;
if force_download || !isfile(data_file)
    println("Downloading data:")
    println("  url  = $data_url")
    println("  file = $data_file")

    open(data_file, "w") do f
        req = HTTP.request(:GET, data_url)
        @assert req.status == 200 "Error while downloading"
        write(f, req.body)
    end

    open("download.stamp", "w") do f
        write(f, string(today()))
    end
else
    println("Using local data:")
    println("  file = $data_file")
    println("  downloaded = ", readline("download.stamp"))
end
```


## Lecture

```julia; echo=false; results="hidden"
@info "* Parsing data"
```

Les données d'entrée sont stockées au format CSV, et contiennent 44 lignes
d'informations préliminaires que nous reproduisons ici, et qui seront sautées
lors de la lecture des données.

```julia; wrap=false
skip = 44
open(data_file) do f
    for _ in 1:skip; println(readline(f)); end
end
```

Le fichier est structuré en deux colonnes :
- `date` : date de la mesure
- `co2` : concentration en CO2 (en ppm molaires)

```julia; results="hidden"
data_raw = CSV.read(data_file; skipto=skip+1, header=[:date, :co2])
```

L'examen des premières et dernières lignes de données révèle qu'elles couvrent
la période de fin mars 1958 jusqu'à nos jours.
```julia; results="raw"
info(data_raw)
```


## Vérification des données manquantes

```julia; echo=false; results="hidden"
@info "* Checking for missing values"
```

Les relevés étant hebdomadaires, l'écart entre deux dates successives du jeu de
données devrait être de 7 jours. Un point manquant provoque un écart de 14
jours, ce qui devrait être rattrappable dans le reste de l'analyse ; au delà, il
faudra se poser des questions sur le traitement à apporter.
```julia
dates = data_raw.date
for i in 2:length(dates)
    if dates[i]-dates[i-1] > 14Days
        println("Missing data: ",
                dates[i-1], " - ", dates[i],
                "  (", dates[i]-dates[i-1], ")")
    end
end
```
Il y a 12 périodes durant lesquelles les données sont manquantes, dont une en
particulier ayant duré 19 semaines en 1964. Le traitement devra en tenir compte.


## Aperçu global des données

```julia; echo=false; results="hidden"
@info "* Plotting raw data"
```

Une visualisation de l'ensemble des données semble montrer une augmentation
tendancielle de la concentration en CO2, à laquelle se superpose une oscillation
à plus haute fréquence.
```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date, data_raw.co2, label="mesures")
disp()
```

Un zoom sur les dernières années laisse penser que le cycle court se produit sur
une période annuelle, avec un minimum local atteint chaque année autour du mois
d'octobre.
```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date[end-150:end], data_raw.co2[end-150:end], label="mesures")
disp()
```

# Analyse

Dans cette analyse, nous allons tenter de séparer ces deux composantes :
composante tendancielle "lisse" et composante oscillante de période
annuelle. Plus précisément, en notant $C$ la concentration en CO2 et $t$ le
temps, nous cherchons une approximation des mesures sous la forme :
$$
    C(t) \approx \theta(t) + \phi(t),
$$
où $\phi(t)$ est une fonction périodique, de période 1 an, donnant la *forme* de
la variation de la concentration en CO2 sur des échelles de temps courtes. $\theta(t)$ est
une fonction très régulière, idéalement un polynôme d'ordre bas, donnant la
*tendance* de la variation de concentration en CO2 en temps long.

La démarche que nous suivons est globalement la suivante : le jeu de données va
être découpé en périodes annuelles. Chacune de ces périodes annuelles sera
traitée (indépendamment des autres) afin d'en extraire une composante lisse, et
une composante périodique.

Si notre hypothèse est correcte, les composantes périodiques de chaque année
devraient être relativement comparables les unes aux autres, et pouvoir être
approchées par leur moyenne : ce comportement annuel moyen constituera notre
composante oscillante $\phi(t)$.

On pourra ensuite obtenir la composante lisse tendancielle en éliminant la
composante oscillatoire moyenne du signal d'origine :
$$
   \theta(t) = C(t)-\phi(t).
$$

La périodicité des données (hebdomadaire) n'étant que peu adaptée à un découpage
annuel, nous allons commencer par interpoler les données à une maille
journalière. Ceci nous permettra de découper le jeu de données en années.


## Travaux sur les dates

Il est plus simple d'interpoler entre deux nombres qu'entre deux dates. Dans la
suite, nous adopterons une convention selon laquelle chaque date peut être
représentée par le nombre de jours qui la sépare de la première mesure :
```julia; results="hidden"
date2num(d::Date) = Dates.value(d - data_raw.date[1])
num2date(n::Int)  = data_raw.date[1] + n * Days
```
Par exemple, pour les premières mesures :
```julia; echo=false
DataFrame(date=data_raw.date[1:3], date_num=date2num.(data_raw.date[1:3]))
```

Afin de comparer des données année par année, nous allons aussi enrichir les
données avec de nouvelles représentations de la date : une date peut être
décomposée comme un couple `(year, day)` dans lequel
- `year` représente l'année
- `day` représente l'indice du jour dans l'année (entre 0 et 365).

Dans ce formalisme, le 1er janvier 2020 est représenté par le couple
`(year=2020, day=0)`. Le 31 décembre 1983 est représenté par le couple
`(year=1983, day=364)`.

```julia
dayinyear(date::Date) = Dates.value(date - Date(year(date)))
dayinyear(num::Int)   = dayinyear(num2date(num))

let d = Date("1983-12-31")
  year(d), dayinyear(d)
end
```

Enfin, prévenons dès maintenant que tout le code de l'analyse fonctionne en
présence d'années bissextiles, mais rien n'a été fait pour les traiter à part :
l'impact, de l'ordre de 1/365 une année sur 4, a été jugé négligeable a priori.


## Interpolation à la maille journalière

```julia; echo=false; results="hidden"
@info "* Interpolating daily data"
```

On construit un interpolateur linéaire basé sur les mesures de CO2 en fonction
de la "date numérique". C'est le paquet Julia
[`DataInterpolations`](https://github.com/PumasAI/DataInterpolations.jl) qui se
charge d'effectuer le gros du travail.

```julia
interp = DI.LinearInterpolation(data_raw.co2, date2num.(data_raw.date));
```

Nous allons profiter de la construction des données interpolées pour gérer le
problème des données manquantes : nous n'interpolerons aucune donnée dans les
"trous" de 3 semaines ou plus. La construction de ce nouveau jeu de données
interpolées est aussi l'occasion d'enrichir les formats de représentation des
dates. Nous avons maintenant 5 colonnes dans notre jeu de données interpolé :
- `date` : date identifiant le jour de la mesure (ou de la valeur interpolée)
- `co2` : valeur de la mesure de CO2 (ou de l'interpolation)
- `date_num` : date convertie en nombre
- `year` : composante de la date identifiant l'année
- `day` : composante de la date identifiant le jour dans l'année (entre 0 et 365)

```julia; results="raw"
data_interp = DataFrame(date_num=Int[], co2=Float64[], date=Date[], year=Int[], day=Int[])
for i in 2:length(dates)
    if dates[i]-dates[i-1] > 14Days
        # pas d'interpolation : seule la date de gauche est inclue
        range = dates[i-1] => dates[i-1]
    else
        # Interpolation entre le début de la période et
        # la fin (exclue car traitée en tant que début de
        # la prochaine période)
        range = dates[i-1] => dates[i]- 1Days
    end

    # Pour chaque jour dans la période considérée,
    # on ajoute une nouvelle ligne de données en interpolant
    for date_num in date2num(first(range)):date2num(last(range))
        date = num2date(date_num)
        push!(data_interp, (date_num = date_num,
                            date = date,
                            year = year(date),
                            day  = dayinyear(date),
                            co2  = interp(date_num)))
    end
end

info(data_interp)
```

En zoomant sur les données interpolées autour de l'une des périodes de données
manquantes, on observe bien le résultat attendu : une interpolation linéaire
journalière lorsque les données sont disponibles, mais aucune interpolation
lorsque les données sont manquantes.

```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]")

idx1 = (data_interp.date .> Date("2006-01-20")) .& (data_interp.date .< Date("2006-03-15"))
plot!(data_interp.date[idx1], data_interp.co2[idx1], seriestype=:scatter,
      label="valeurs interpolées")

idx2 = (data_raw.date .> Date("2006-01-20")) .& (data_raw.date .< Date("2006-03-15"))
plot!(data_raw.date[idx2], data_raw.co2[idx2], seriestype=:scatter,
      label="mesures")
disp()
```

Enfin, nous ne gardons pour notre analyse que des années complètes. D'ailleurs,
afin de tester la validité de nos résultats, nous n'allons réaliser l'analyse
que sur une fraction des données ; nous garderons les 5 dernières années comme
données de test. Les données sur lesquelles portera l'analyse couvrent donc la
période 1959-2014.

```julia; results="raw"
firstyear = minimum(data_interp.year)+1 # Première année incomplète
lastyear  = maximum(data_interp.year)-6 # Dernière année incomplète + 5 années de test
idx = (data_interp.year.>=firstyear) .& (data_interp.year.<=lastyear)
data = data_interp[idx, :];

info(data)
```

Sur ces années complètes, la composante `day` de la date devrait être
équirépartie entre 0 et 365, ce qui est globalement le cas. Les données
manquantes n'ont donc pas d'impact significatif de ce point de vue là.

```julia; echo=false; results="raw"
plot(xlabel="Jour de l'année", ylabel="Nombre de valeurs")
histogram!(data.day, bins=0:31:365, label=nothing)
disp()
```

## Analyse des variations annuelles

```julia; echo=false; results="hidden"
@info "* Fitting yearly oscillations"
```

Pour chaque année, on commence par tenter d'extraire la composante oscillante de
la mesure. Si l'on note $C_a(d)$ la concentration en CO2 le jour numéro $d$ de
l'année $a$, on cherche à écrire :
$$
  \forall d\in{0\ldots365}, \quad
  C_a(d)
  = \underbrace{\theta_a(d)}_{\text{tendance locale}}
  + \underbrace{\phi_a(d)}_{\text{forme locale}}
$$
ce qui correspond à une version locale (pour l'année $a$) de l'expression
globale cherchée pour la variation de concentration en CO2 : $\theta_a$ et
$\phi_a$ donnent respectivement la tendance et la forme de la concentration pour
l'année $a$.

Cherchant une tendance d'ordre aussi bas que possible, nous supposons que
$\theta_a$ peut être approchée par un modèle affine à cette échelle de temps
courts :
$$
  \theta_a(d) = \alpha_a + \beta_a \, d.
$$

Pour que la forme $\phi_a$ soit périodique, il faut que cette fonction prenne la
même valeur en début d'année qu'en fin d'année : $\phi(0)=\phi(365)$. Afin de
définir la constante $\alpha$ de manière unique, on fixe de plus :
$$
  \phi(0)=\phi(365)=0.
$$

On obtient donc
$$
\begin{align*}
  \alpha_a  &= C_a(0),                      \\[1em]
  \beta_a   &= \frac{C_a(365)-C_a(0)}{365},   \\[1em]
  \phi_a(d) &= C_a(d) - \alpha_a - \beta_a d.
\end{align*}
$$

En pratique, plutôt que des valeurs ponctuelles $C_a(0)$ et $C_a(365)$, on prend
plutôt des valeurs (notées $C_0$ et $C_1$ dans le code) moyennées sur les 7
premiers et 7 derniers jours de l'année.

```julia
data.phi   = copy(data.co2)
data.alpha = zero(data.co2)
data.beta  = zero(data.co2)
for year in unique(data.year)
    idx = (data.year .== year)
    C₀ = data.co2[idx .& (data.day .< 7)  ] |> mean
    C₁ = data.co2[idx .& (data.day .> 358)] |> mean
    α = C₀
    β = (C₁ - C₀) / 365

    data.alpha[idx] .= α
    data.beta[idx]  .= β
    data.phi[idx]   .-= α .+ β*data.day[idx]
end
```

Examinons par exemple l'effet de ce traitement sur les données interpolées de
l'année 1982. On voit, sur la figure du haut, les mesures brutes ainsi que la
tendance locale (affine). Sur la figure du bas, la composante périodique locale
vérifie bien les contraintes demandées, avec ses valeurs nulles aux bords.

```julia; echo=false; results="raw"
let df = data[data.year.==1982, :]
   p1 = plot(ylabel="CO2 [ppm]")
   plot!(p1, df.date, df.co2,
         label="mesures", linewidth=2)
   plot!(p1, df.date, df.alpha.+df.beta.*df.day,
         label="alpha + beta * day", linewidth=2)


   p2 = plot(xlabel="date", ylabel="CO2 [ppm]")
   plot!(p2, df.date, df.phi,
         label="phi", linecolor="green", linewidth=2)

   l = @layout [a ; b]
   plot(p1, p2, layout=l)
   disp()
end
```

Pour que notre décomposition soit valide, il faut que les composantes
périodiques locales de chaque année se ressemblent, au point de pouvoir être
représentées par leur moyenne. Nous calculons donc cette moyenne (et
l'écart-type associé) pour toutes les années du jeu de données d'étude.

```julia; results="raw"
avg = by(data, :day, :phi=>mean, :phi=>std)
info(avg)
```

Et nous traçons l'ensemble des composantes oscillantes locales aux côtés de cette moyenne.

```julia; echo=false; results="raw"
plot(xlabel="date", ylabel="CO2 [ppm]")
plot!(data.day, data.phi, seriestype=:scatter, label="phi")
plot!(avg.day, avg.phi_mean, linewidth=3, label="phi (moyenne)")
plot!(avg.day, avg.phi_mean.-2avg.phi_std, linecolor="red", linewidth=2, label="phi (IC 95%)")
plot!(avg.day, avg.phi_mean.+2avg.phi_std, linecolor="red", linewidth=2, label=nothing)
plot!(avg.day, avg.phi_std, linewidth=2, label="phi (écart-type)")
disp()
```

Comme on pouvait s'y attendre, il reste une forte variabilité d'année en
année. Notons toutefois que l'écart-type, bien qu'important, reste un ordre de
grandeur en dessous des valeurs extrêmes, ce qui permet d'espérer que ce profil
moyen est suffisamment représentatif pour permettre d'extraire la composante
tendancielle lisse des données.


## Analyse des variations tendancielles

```julia; echo=false; output="hidden"
@info "* Fitting underlying trend"
```

Nous sommes maintenant prêts à extraire la composante tendancielle des
mesures. Il suffit pour cela de retrancher aux données brutes la composante
oscillante moyenne : d'après notre modèle, on a en effet
$$
  \theta(t) = C(t) - \phi(t).
$$

```julia
data = join(data, avg, on=:day)
data.theta = data.co2 .- data.phi_mean;
```

Même s'il reste des oscillations locales, nous constatons tout de même que la
composante tendancielle est devenue suffisamment lisse pour récupérer une forme
de monotonie.

```julia; echo=false; results="raw"
plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data.date, data.co2,   label="mesures")
plot!(data.date, data.theta, label="theta")
disp()
```

Nous allons maintenant tenter de caractériser la tendance sous-jacente. Au vu de
la courbe (convexe), nous proposons un modèle quadratique de la forme :
$$
\theta(t) = \alpha + \beta t + \gamma t^2
$$

Le paquet Julia [`GLM`](https://github.com/JuliaStats/GLM.jl) permet de *fitter*
des modèles linéaires (ou modèles linéaires généralisés). Nous l'utilisons ici
pour estimer les paramètres de notre modèle. Afin de ne pas considérer un nombre
artificiellement élevé de points de mesure, l'estimation est réalisée sur un
sous-échantillonnage du jeu de données, comprenant un point (moyen) par an, ce
qui correspond à l'échelle de temps que nous n'avons pas encore capturé dans la
composante oscillante.

```julia; wrap=false
sample = by(data, :year, theta = :theta=>mean, date_num = :date_num=>mean)
model  = GLM.lm(@formula(theta ~ date_num + date_num^2), sample)
```

Le modèle obtenu semble correspondre assez bien aux données, avec une
incertitude de l'ordre de quelques pourcents sur l'estimation des paramètres
$\alpha$ et $\beta$.

```julia; results="hidden"
# Paramètres du modèle et intervalle de confiance à 95%
α,  β,  γ  = coef(model)
α1, β1, γ1 = coef(model)-2*stderror(model)
α2, β2, γ2 = coef(model)+2*stderror(model)

# Taux d'accroissement
τ(d) = (β + γ * date2num(d)) * 365

τ1 = τ(Date("1959"))
τ2 = τ(Date("2015"))
```

Ce modèle caractérise en particulier une hausse tendancielle de la concentration
de CO2 atmosphérique de l'ordre de $\tau_1$ = `j @printf "%.2f" τ1` ppm/an en
1959 et $\tau_2$ = `j @printf "%.2f" τ2` ppm/an en 2015. Si la valeur de 1959
correspond assez bien au taux de 0.75 ppm/an déterminé dans `j monroe2015()`, la
valeur que nous trouvons pour 2015 est sous-évaluée de plus de 30% par rapport
aux 2.25 ppm/an calculés par Monroe.

Par ailleurs, s'il est clair que la tendance est à l'augmentation, on voit
toutefois que l'incertitude sur les paramètres n'est pas complètement
négligeable. L'incertitude sur $\beta$ est en particulier de nature à engendrer
une perte de prédictibilité du modèle en temps long.

```julia; echo=false; results="raw"
data.theta_pred = predict(model, data)

plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date, data_raw.co2, label="mesures")
plot!(data.date, data.theta_pred, label="modèle theta", linewidth=3)
plot!(data.date, α .+ β1 * data.date_num .+ γ*data.date_num.^2,
      linecolor="red", label="modèle theta (IC 95% sur beta)")
plot!(data.date, α .+ β2 * data.date_num .+ γ*data.date_num.^2,
      linecolor="red", label=nothing)
disp()
```

## Reconstruction du signal complet et prédiction

```julia; echo=false; output="hidden"
@info "* Reconstructing and predicting complete signal"
```

Nous avons maintenant tous les éléments nécessaires afin de reconstruire
l'évolution de la concentration en CO2 selon notre modèle :
$$
  C(t) = \theta(t) + \phi(t).
$$

Les 5 dernières années du jeu de données brut n'ont pas servi à calibrer notre
modèle ; nous allons les utiliser comme données de test. Nous extrapolons aussi
les valeurs des 5 prochaines années. Pour toute cette période, nous
reconstruisons les valeurs de concentration en CO2 tous les 10 jours.

```julia; results="raw"
date_num = date2num(Date(lastyear)):10:date2num(today()+5Years)
prediction = DataFrame(date_num=date_num,
                       date=num2date.(date_num),
                       day=dayinyear.(date_num))
info(prediction)
```

Le modèle `GLM` précédemment calibré est utilisé pour prédire $\theta(t)$. On y
ajoute la fonction de forme annuelle $\phi(t)$ par périodicité.

```julia; results="raw"
prediction.theta = predict(model, prediction)
prediction = join(prediction, avg, on=:day)

prediction.co2 = prediction.theta .+ prediction.phi_mean
info(prediction)
```

On voit que la forme annuelle semble bien reproduite sur les 5 premières années,
pour lesquelles il est possible de comparer les prédictions avec les mesures
réelles. En revanche, la tendance ne colle que sur les deux premières années
d'extrapolation ; on observe un décalage significatif et croissant ensuite.

```julia; echo=false; results="raw"
idx = data_raw.date .> Date(lastyear)
plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date[idx], data_raw.co2[idx], label="mesures")
plot!(prediction.date, prediction.theta,     label="tendance", linewidth=2)
plot!(prediction.date, prediction.co2,       label="prédiction", linewidth=2)
disp()
```

## Pour aller plus loin

```julia; echo=false
@info "* Studying alternate model"
```

En utilisant des techniques de calage de modèle plus complexes, il est possible
de mieux caractériser la tendance en temps longs :

```julia; wrap=false
model2 = glm(@formula(theta ~ date_num + date_num^2), sample,
             InverseGaussian(), InverseSquareLink())
```

En reprenant l'analyse précédente, ce nouveau modèle donne les prédictions
suivantes, qui collent quasi-parfaitement aux mesures dans la période de test :

```julia; echo=false; results="raw"
date_num = date2num(Date(lastyear)):10:date2num(today()+5Years)
prediction2 = DataFrame(date_num=date_num,
                        date=num2date.(date_num),
                        day=dayinyear.(date_num))

prediction2.theta = predict(model2, prediction2)
prediction2 = join(prediction2, avg, on=:day)
prediction2.co2 = prediction2.theta .+ prediction2.phi_mean

idx = data_raw.date .> Date(lastyear)
plot(xlabel="Date", ylabel="CO2 [ppm]")
plot!(data_raw.date[idx], data_raw.co2[idx], label="mesures")
plot!(prediction2.date, prediction2.theta,   label="tendance", linewidth=2)
plot!(prediction2.date, prediction2.co2,     label="prédiction", linewidth=2)
disp()
```

Je n'ai aucune idée de la signification qu'il est possible de donner à ces
valeurs de paramètres, aussi est-il sans doute préférable de conserver, au moins
dans un premier temps, l'interprétation donnée par le modèle quadratique simple
(quoi que donnant des prédictions plus éloignées des données).


# Conclusions

Dans cette étude, nous avons tenté de proposer une analyse (reproductible) de la
courbe de Keeling, qui suit les variations de la concentration atmosphérique en
CO2 depuis 1959 à Hawaï. La courbe est décomposée en une composante tendancielle
à laquelle se superposent des variations périodiques annuelles. La coomposante
tendancielle montre clairement un comportement d'augmentation sur le long
terme. Une caractérisation de cette tendance sous forme de modèle quadratique
fait apparaître des taux d'accroissement de l'ordre de 0.8 ppm/an en 1959, qui
s'accélèrent pour monter à environ 1.5 ppm/an en 2015. Les taux d'accroissement
trouvés ici peuvent être comparés aux taux de 0.75 ppm/an (1959) et 2.25 ppm/an
(2015) trouvés dans la littérature `j monroe2015()`. Il serait intéressant de
pousser la comparaison entre ces deux travaux plus avant, afin de comprendre
l'origine des écarts significatifs (de l'ordre de 30%) sur les taux
d'accroissements de 2015.

Par ailleurs, nous avons dans cette étude proposé une validation du modèle
utilisant les 5 dernières années de mesures comme données de test. Une
extrapolation sur les 5 prochaines années (2020-2025) est aussi proposée. Il
convient de noter ici que l'utilisation de techniques poussées de calage de
paramètre semble de prime abord donner des estimations plus prédictives. Une
autre perspective de ce travail pourrait être d'interpréter et valider
l'utilisation de telles techniques dans ce contexte.

<!-- Local Variables: -->
<!-- mode: markdown -->
<!-- ispell-local-dictionary: "french" -->
<!-- End: -->