From 275baf11d80ca8d1db4de35667bd02ac6d4761bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20F=C3=A9votte?= Date: Sat, 11 Apr 2020 01:35:34 +0200 Subject: [PATCH] Ex 03-3 (suite) --- module3/exo3/exercice.html | 243 +++++++++++++++++--------------- module3/exo3/exercice.jmd | 277 +++++++++++++++++++++++++------------ module3/exo3/utils.jl | 61 ++++++++ 3 files changed, 378 insertions(+), 203 deletions(-) create mode 100644 module3/exo3/utils.jl diff --git a/module3/exo3/exercice.html b/module3/exo3/exercice.html index 28dd2b4..8ede906 100644 --- a/module3/exo3/exercice.html +++ b/module3/exo3/exercice.html @@ -662,7 +662,7 @@ div.title {text-align: center;}

Gestion des dépendances

-

Environnement

+

Environnement

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

@@ -700,7 +700,7 @@ Environment: -

Chargement des dépendances

+

Chargement des dépendances

Tant que nous y sommes, profitons en pour charger dès maintenant les paquets dont nous aurons besoin.

@@ -720,13 +720,14 @@ Environment: +

Données d'entrée

Nos données d'entrée proviennent du programme Scripps CO2. Nous fondons l'analyse sur le jeu de données contenant des observations hebdomadaires.

-

Téléchargement

+

Téléchargement

-

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 version des données utilisées pour l'analyse, et stockée dans git aux côtés du présent document.

-

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

+

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.

@@ -768,7 +769,7 @@ Using local data:
-

Lecture

+

Lecture

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.

@@ -848,29 +849,16 @@ Using local data:

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.

- 
 info(data_raw)
-
-3156×2 DataFrame
-│ Row  │ date       │ val     │
-│      │ Date       │ Float64 │
-├──────┼────────────┼─────────┤
-│ 1    │ 1958-03-29 │ 316.19  │
-│ 2    │ 1958-04-05 │ 317.31  │
-│ 3    │ 1958-04-12 │ 317.69  │
-⋮
-│ 3153 │ 2020-01-04 │ 413.19  │
-│ 3154 │ 2020-01-11 │ 413.39  │
-│ 3155 │ 2020-01-25 │ 413.36  │
-│ 3156 │ 2020-02-01 │ 413.99  │
-
+3156 rows × 2 columns +
dateval
DateFloat64
11958-03-29316.190000
21958-04-05317.310000
31958-04-12317.690000
...
31542020-01-11413.390000
31552020-01-25413.360000
31562020-02-01413.990000
-

Vérification des données manquantes

+

Vérification des données manquantes

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.

@@ -879,7 +867,7 @@ Using local data: 
 dates = data_raw.date
 for i in 2:length(dates)
-    if dates[i]-dates[i-1] > Dates.Day(14)
+    if dates[i]-dates[i-1] > 14Days
         println("Missing data: ",
                 dates[i-1], " - ", dates[i],
                 "  (", dates[i]-dates[i-1], ")")
@@ -905,29 +893,37 @@ Missing data: 2012-09-29 - 2012-10-20  (21 days)
 
 
 
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

+
Aperçu global des données

 
 
 
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.

 
-
+
 
 
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.

 
-
+
 
 
Analyse

-
Dans cette analyse, nous allons tenter de séparer ces deux composantes : composante tendancielle "lisse" et composante oscillante de période annuelle.

+
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. On peut ensuite obtenir la composante lisse tendancielle en éliminant la composante oscillatoire moyenne du signal d'origine.

+
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

+
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 :

 
 
 
 date2num(d::Date) = Dates.value(d - data_raw.date[1])
-num2date(n::Int)  = data_raw.date[1] + Day(n)
+num2date(n::Int)  = data_raw.date[1] + n * Days
 

 
 
@@ -963,7 +959,7 @@ Missing data: 2012-09-29 - 2012-10-20  (21 days)
 
 
 
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

+
Interpolation à la maille journalière

 
 
 
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 qui se charge d'effectuer le gros du travail.

@@ -993,14 +989,14 @@ Missing data: 2012-09-29 - 2012-10-20  (21 days)
 
 data_interp = DataFrame(date_num=Int[], val=Float64[], date=Date[], year=Int[], day=Int[])
 for i in 2:length(dates)
-    if dates[i]-dates[i-1] > Dates.Day(14)
+    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]-Day(1)
+        range = dates[i-1] => dates[i]- 1Days
     end
 
     # Pour chaque jour dans la période considérée,
@@ -1015,78 +1011,67 @@ Missing data: 2012-09-29 - 2012-10-20  (21 days)
     end
 end
 
-# 3 dernières lignes, pour vérification
-last(data_interp, 3)
+info(data_interp)
 

 
 
+22132 rows × 5 columns
+
date_numvaldateyearday
Int64Float64DateInt64Int64
10316.1900001958-03-29195887
21316.3500001958-03-30195888
32316.5100001958-03-31195889
...
2213022586413.7200002020-01-29202028
2213122587413.8100002020-01-30202029
2213222588413.9000002020-01-31202030

 
-
3 rows × 5 columns
date_numvaldateyearday
Int64Float64DateInt64Int64
122586413.722020-01-29202028
222587413.812020-01-30202029
322588413.92020-01-31202030

 
 
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.

 
-
+
 
-
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.

+
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.

 
 
 
 firstyear = minimum(data_interp.year)+1
 lastyear  = maximum(data_interp.year)-6
 idx = (data_interp.year.>=firstyear) .& (data_interp.year.<=lastyear)
-data = data_interp[idx, :];
-

-
-
+data = data_interp[idx, :];
 
-
Les données sur lesquelles portera l'analyse couvrent donc la période 1959-2014.

-
-
-
-info(data)
+info(data)
 

 
 
-
-20100×5 DataFrame
-│ Row   │ date_num │ val     │ date       │ year  │ day   │
-│       │ Int64    │ Float64 │ Date       │ Int64 │ Int64 │
-├───────┼──────────┼─────────┼────────────┼───────┼───────┤
-│ 1     │ 278      │ 315.231 │ 1959-01-01 │ 1959  │ 0     │
-│ 2     │ 279      │ 315.236 │ 1959-01-02 │ 1959  │ 1     │
-│ 3     │ 280      │ 315.24  │ 1959-01-03 │ 1959  │ 2     │
-⋮
-│ 20097 │ 20728    │ 399.14  │ 2014-12-28 │ 2014  │ 361   │
-│ 20098 │ 20729    │ 399.27  │ 2014-12-29 │ 2014  │ 362   │
-│ 20099 │ 20730    │ 399.4   │ 2014-12-30 │ 2014  │ 363   │
-│ 20100 │ 20731    │ 399.53  │ 2014-12-31 │ 2014  │ 364   │
-

+20100 rows × 5 columns
+
date_numvaldateyearday
Int64Float64DateInt64Int64
1278315.2314291959-01-0119590
2279315.2357141959-01-0219591
3280315.2400001959-01-0319592
...
2009820729399.2700002014-12-292014362
2009920730399.4000002014-12-302014363
2010020731399.5300002014-12-312014364

 
 
-
Sur ces années complètes, la composante day de la date devrait être globalement é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à.

+
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à.

 
-
+
 
-
Analyse des variations annuelles

+
Analyse des variations annuelles

 
 
-
Pour chaque année, on commence par tenter d'extraire la composante oscillante de la mesure. En supposant que cette composante est périodique, il faut que la valeur de cette composante en début d'année soit égale à celle de fin d'année. On va aussi supposer que la composante tendancielle varie linéairement en cours d'année. Si on note $C$ la concentration en CO2 et $d$ le jour, on cherche à écrire :

+
Pour chaque année, on commence par tenter d'extraire la composante oscillante de la mesure. Si l'on note $C_a$ la concentration en CO2 durant l'année $a$ et $d$ un jour de cette année, on cherche à écrire :

 
\[
   \forall d\in{0\ldots365}, \quad
-  C(d)
-  = \underbrace{\alpha + \beta d}_{\text{tendance locale}}
-  + \underbrace{\phi(d)}_{\text{composante oscillante}}
+  C_a(d)
+  = \underbrace{\theta_a(d)}_{\text{tendance locale}}
+  + \underbrace{\phi_a(d)}_{\text{forme locale}}
+\]

+
ce qui correspond à une version locale 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.
 \]

-
avec $\phi(0)=\phi(365)$ pour garantir la périodicité de la composante oscillante. 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  &= C(0),                      \\[1em]
-  \beta   &= \frac{C(365)-C(0)}{365},   \\[1em]
-  \phi(d) &= C(d) - \alpha - \beta d.
+  \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(0)$ et $C(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.

+
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.

 
 
 
@@ -1110,55 +1095,63 @@ Missing data: 2012-09-29 - 2012-10-20  (21 days)
 
 
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.

 
-
+
 
 
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 pour toutes les années du jeu de données d'étude.

 
 
 
-avg = by(data, :day, :phi=>mean, :phi=>std);
+avg = by(data, :day, :phi=>mean, :phi=>std)
+info(avg)
 

 
 
+366 rows × 3 columns
+
dayphi_meanphi_std
Int64Float64Float64
10-0.1022120.125698
21-0.0738840.088030
32-0.0439490.048068
...
3643630.0589410.055141
3653640.0915550.077335
3663650.1266760.100662

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

 
-
+
 
 
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

+
Analyse des variations tendancielles

 
 
-
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.

+
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).
+\]

 
 
 
 data = join(data, avg, on=:day)
-data.smooth = data.val .- data.phi_mean;
+data.theta = data.val .- 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.

 
-
+
 
-
Nous allons maintenant tenter de caractériser la tendance sous-jacente. Au vu de la courbe (concave), nous proposons un modèle quadratique de la forme :

+
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 :

 
\[
-C(t) = \alpha + \beta t + \gamma t^2
+\theta(t) = \alpha + \beta t + \gamma t^2
 \]

+
Le paquet Julia GLM 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. L'estimation est réalisée sur un sous-échantillonnage du jeu de données, comprenant un point (moyen) par an.

 
 
 
-model = lm(@formula(smooth ~ date_num + date_num^2),
-            by(data, :year, smooth = :smooth=>mean, date_num = :date_num=>mean))
+sample = by(data, :year, theta = :theta=>mean, date_num = :date_num=>mean)
+model  = GLM.lm(@formula(theta ~ date_num + date_num^2), sample)
 

 
 
 
 StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}
 
-smooth ~ 1 + date_num + :(date_num ^ 2)
+theta ~ 1 + date_num + :(date_num ^ 2)
 
 Coefficients:
 ───────────────────────────────────────────────────────────────────────────────────────
@@ -1171,7 +1164,7 @@ date_num ^ 2    9.10059e-8  2.66606e-9    34.135     <1e-37    8.56585e-8
 

 
 
-
Ce genre de modèle semble correspondre assez bien aux données, avec un assez bon niveau de confiance sur l'estimation des paramètres $\alpha$ et $\beta$.

+
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$.

 
 
 
@@ -1183,62 +1176,86 @@ date_num ^ 2    9.10059e-8  2.66606e-9    34.135     <1e-37    8.56585e-8
 
 
 
En première approximation, on a en particulier une hausse tendancielle de la concentration de CO2 atmosphérique de l'ordre de $\beta$ = 2.22e-03 ppm/jour ($\pm$5.21%), qui se traduit en une augmentation annuelle comprise entre 0.77 et 0.85 ppm/an avec un niveau de confiance de 95%.

-
On voit toutefois que l'incertitude sur $\beta$ en particulier est de nature à engendrer des

+
On voit toutefois que l'incertitude sur les paramètres (en particulier $\beta$) est de nature à engendrer une perte de prédictibilité du modèle en temps long.

+
+
+
+
Reconstruction du signal complet et prédiction

 
-
+
+
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.

 
 
 
-date_num = date2num(Date(lastyear)):10:date2num(today()+Year(5))
+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)
+

 
-prediction.smooth = predict(model, prediction)
-prediction = join(prediction, avg, on=:day)
 
-prediction.val = prediction.smooth .+ prediction.phi_mean
+412 rows × 3 columns
+
date_numdateday
Int64DateInt64
1203672014-01-010
2203772014-01-1110
3203872014-01-2120
...
410244572025-03-1472
411244672025-03-2482
412244772025-04-0392

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

+
 
-idx = data_raw.date .> Date(lastyear)
-plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
-plot!(data_raw.date[idx], data_raw.val[idx], label="measurements")
-plot!(prediction.date, prediction.smooth,    label="smooth model")
-plot!(prediction.date, prediction.val,       label="predicted value")
+
+prediction.theta = predict(model, prediction)
+prediction = join(prediction, avg, on=:day)
+
+prediction.val = prediction.theta .+ prediction.phi_mean
+info(prediction)
 

 
 
-
+412 rows × 7 columns
+
date_numdatedaythetaphi_meanphi_stdval
Int64DateInt64Float64?Float64Float64Float64
1203672014-01-010396.753169-0.1022120.125698396.650956
2203772014-01-1110396.8124120.1981870.220453397.010599
3203872014-01-2120396.8716730.4034950.374229397.275168
...
410244572025-03-1472422.5021751.6335060.561281424.135681
411244672025-03-2482422.5688622.0232960.459201424.592158
412244772025-04-0392422.6355682.5383720.467285425.173939

+
+
+
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.

 
-
Pour aller plus loin

+
+
+
Pour aller plus loin

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

 
 
 
-model = glm(@formula(smooth ~ date_num),
-            by(data, :year, smooth = :smooth=>mean, date_num = :date_num=>mean),
-            InverseGaussian(), InverseSquareLink())
+model2 = glm(@formula(theta ~ date_num + date_num^2), sample,
+             InverseGaussian(), InverseSquareLink())
 

 
 
 
 StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},InverseGaussian{Float64},InverseSquareLink},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}
 
-smooth ~ 1 + date_num
+theta ~ 1 + date_num + :(date_num ^ 2)
 
 Coefficients:
-─────────────────────────────────────────────────────────────────────────────────────
-                 Estimate  Std. Error   z value  Pr(>|z|)     Lower 95%     Upper 95%
-─────────────────────────────────────────────────────────────────────────────────────
-(Intercept)   1.02817e-5   1.31515e-8   781.789    <1e-99   1.02559e-5    1.03075e-5
-date_num     -1.91994e-10  1.0025e-12  -191.515    <1e-99  -1.93959e-10  -1.90029e-10
-─────────────────────────────────────────────────────────────────────────────────────
+────────────────────────────────────────────────────────────────────────────────────────
+                  Estimate   Std. Error    z value  Pr(>|z|)     Lower 95%     Upper 95%
+────────────────────────────────────────────────────────────────────────────────────────
+(Intercept)    1.02097e-5   1.68812e-8   604.796      <1e-99   1.01766e-5    1.02428e-5
+date_num      -1.73581e-10  3.46014e-12  -50.166      <1e-99  -1.80363e-10  -1.668e-10
+date_num ^ 2  -8.3699e-16   1.52956e-16   -5.47208    <1e-7   -1.13678e-15  -5.37201e-16
+────────────────────────────────────────────────────────────────────────────────────────
 

 
 
+
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 :

 
-
-closeall()
-

+
 
+
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).

 
 
 
@@ -1257,7 +1274,7 @@ date_num     -1.91994e-10  1.0025e-12  -191.515    <1e-99  -1.93959e-10  -1.9
           

           Published from exercice.jmd using
           Weave.jl
-           on 2020-04-10.
+           on 2020-04-11.
           

 
 
diff --git a/module3/exo3/exercice.jmd b/module3/exo3/exercice.jmd
index bcfb6c0..3baaee4 100644
--- a/module3/exo3/exercice.jmd
+++ b/module3/exo3/exercice.jmd
@@ -7,7 +7,7 @@ date: avril 2020
 # Gestion des dépendances
 
 
-## Environnement
+### Environnement
 
 ```julia; echo=false; results="hidden"
 @info "Instantiating project"
@@ -37,7 +37,7 @@ Pkg.instantiate()
 ```
 
 
-## Chargement des dépendances
+### Chargement des dépendances
 
 ```julia; echo=false; results="hidden"
 @info "Loading dependencies"
@@ -57,6 +57,9 @@ using Statistics
 using Plots; gr()
 ```
 
+```julia; echo=false; results="hidden"
+include("utils.jl")
+```
 
 # Données d'entrée
 
@@ -65,18 +68,20 @@ CO2](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.htm
 fondons l'analyse sur le jeu de données contenant des observations hebdomadaires.
 
 
-## Téléchargement
+### 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 version des données
-utilisées pour l'analyse, et stockée dans git aux côtés du présent document.
+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 possible de forcer le téléchargement en positionnant la variable
-`force_download=true` afin d'actualiser le jeu de données.
+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"
@@ -107,7 +112,7 @@ end
 ```
 
 
-## Lecture
+### Lecture
 
 ```julia; echo=false; results="hidden"
 @info "Parsing data"
@@ -134,15 +139,12 @@ data_raw = CSV.read(data_file; skipto=skip+1, header=[:date, :val])
 
 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; echo=false
-info(df::DataFrame) = show(IOContext(stdout, :limit=>true, :displaysize=>(15, 100)), df);
-```
-```julia
+```julia; results="raw"
 info(data_raw)
 ```
 
 
-## Vérification des données manquantes
+### Vérification des données manquantes
 
 ```julia; echo=false; results="hidden"
 @info "Checking for missing values"
@@ -155,7 +157,7 @@ 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] > Dates.Day(14)
+    if dates[i]-dates[i-1] > 14Days
         println("Missing data: ",
                 dates[i-1], " - ", dates[i],
                 "  (", dates[i]-dates[i-1], ")")
@@ -166,7 +168,7 @@ 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
+### Aperçu global des données
 
 ```julia; echo=false; results="hidden"
 @info "Plotting raw data"
@@ -177,7 +179,7 @@ tendancielle de la concentration en CO2, à laquelle se superpose une oscillatio
 à plus haute fréquence.
 ```julia; echo=false
 plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
-plot!(data_raw.date, data_raw.val, label="measurements")
+plot!(data_raw.date, data_raw.val, label="mesures")
 ```
 
 Un zoom sur les dernières années laisse penser que le cycle court se produit sur
@@ -185,13 +187,22 @@ une période annuelle, avec un minimum local atteint chaque année autour du moi
 d'octobre.
 ```julia; echo=false
 plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
-plot!(data_raw.date[end-150:end], data_raw.val[end-150:end], label="measurements")
+plot!(data_raw.date[end-150:end], data_raw.val[end-150:end], label="mesures")
 ```
 
 # Analyse
 
 Dans cette analyse, nous allons tenter de séparer ces deux composantes :
-composante tendancielle "lisse" et composante oscillante de période annuelle.
+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
@@ -200,22 +211,28 @@ 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. On peut ensuite obtenir la composante lisse
-tendancielle en éliminant la composante oscillatoire moyenne du signal d'origine.
+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
+### 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] + Day(n)
+num2date(n::Int)  = data_raw.date[1] + n * Days
 ```
 Par exemple, pour les premières mesures :
 ```julia; echo=false
@@ -246,7 +263,7 @@ présence d'années bissextiles, mais rien n'a été fait pour les traiter à pa
 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
+### Interpolation à la maille journalière
 
 ```julia; echo=false; results="hidden"
 @info "Interpolating daily data"
@@ -272,17 +289,17 @@ dates. Nous avons maintenant 5 colonnes dans notre jeu de données interpolé :
 - `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
+```julia; results="raw"
 data_interp = DataFrame(date_num=Int[], val=Float64[], date=Date[], year=Int[], day=Int[])
 for i in 2:length(dates)
-    if dates[i]-dates[i-1] > Dates.Day(14)
+    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]-Day(1)
+        range = dates[i-1] => dates[i]- 1Days
     end
 
     # Pour chaque jour dans la période considérée,
@@ -297,8 +314,7 @@ for i in 2:length(dates)
     end
 end
 
-# 3 dernières lignes, pour vérification
-last(data_interp, 3)
+info(data_interp)
 ```
 
 En zoomant sur les données interpolées autour de l'une des périodes de données
@@ -311,69 +327,80 @@ plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:topleft)
 
 idx1 = (data_interp.date .> Date("2006-01-20")) .& (data_interp.date .< Date("2006-03-15"))
 plot!(data_interp.date[idx1], data_interp.val[idx1], seriestype=:scatter,
-      label="interpolated values")
+      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.val[idx2], seriestype=:scatter,
-      label="measurements")
+      label="mesures")
 ```
 
 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.
+données de test. Les données sur lesquelles portera l'analyse couvrent donc la
+période 1959-2014.
 
-```julia
+```julia; results="raw"
 firstyear = minimum(data_interp.year)+1
 lastyear  = maximum(data_interp.year)-6
 idx = (data_interp.year.>=firstyear) .& (data_interp.year.<=lastyear)
 data = data_interp[idx, :];
-```
 
-Les données sur lesquelles portera l'analyse couvrent donc la période 1959-2014.
-```julia
 info(data)
 ```
 
 Sur ces années complètes, la composante `day` de la date devrait être
-globalement é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à.
+é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
-histogram(data.day, bins=0:31:365, label=nothing)
+plot(xlabel="Jour de l'année", ylabel="Nombre de valeurs")
+histogram!(data.day, bins=0:31:365, label=nothing)
 ```
 
-## Analyse des variations annuelles
+### 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. En supposant que cette composante est périodique, il faut que la
-valeur de cette composante en début d'année soit égale à celle de fin
-d'année. On va aussi supposer que la composante tendancielle varie linéairement
-en cours d'année. Si on note $C$ la concentration en CO2 et $d$ le jour, on
-cherche à écrire :
+la mesure. Si l'on note $C_a$ la concentration en CO2 durant l'année $a$ et
+$d$ un jour de cette année, on cherche à écrire :
 $$
   \forall d\in{0\ldots365}, \quad
-  C(d)
-  = \underbrace{\alpha + \beta d}_{\text{tendance locale}}
-  + \underbrace{\phi(d)}_{\text{composante oscillante}}
+  C_a(d)
+  = \underbrace{\theta_a(d)}_{\text{tendance locale}}
+  + \underbrace{\phi_a(d)}_{\text{forme locale}}
+$$
+ce qui correspond à une version locale 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.
 $$
-avec $\phi(0)=\phi(365)$ pour garantir la périodicité de la composante
-oscillante. 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  &= C(0),                      \\[1em]
-  \beta   &= \frac{C(365)-C(0)}{365},   \\[1em]
-  \phi(d) &= C(d) - \alpha - \beta d.
+  \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(0)$ et $C(365)$, on prend
+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.
 
@@ -402,7 +429,7 @@ vérifie bien les contraintes demandées, avec ses valeurs nulles aux bords.
 ```julia; echo=false
 let df = data[data.year.==1982, :]
    p1 = plot(ylabel="CO2 [ppm]")
-   plot!(p1, df.date, df.val, label="measurements")
+   plot!(p1, df.date, df.val, label="mesures")
    plot!(p1, df.date, df.alpha.+df.beta.*df.day, label="alpha + beta * day")
 
 
@@ -416,10 +443,12 @@ 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 pour toutes les années du jeu de données d'étude.
+représentées par leur moyenne. Nous calculons donc cette moyenne pour toutes les
+années du jeu de données d'étude.
 
-```julia
-avg = by(data, :day, :phi=>mean, :phi=>std);
+```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.
@@ -428,7 +457,9 @@ Et nous traçons l'ensemble des composantes oscillantes locales aux côtés de c
 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_std,  linewidth=1, label="phi (écart-type)")
+plot!(avg.day, avg.phi_mean.-avg.phi_std, linecolor="red", linewidth=1, label="phi (IC 67%)")
+plot!(avg.day, avg.phi_mean.+avg.phi_std, linecolor="red", linewidth=1, label=nothing)
+plot!(avg.day, avg.phi_std, linewidth=1, label="phi (écart-type)")
 ```
 
 Comme on pouvait s'y attendre, il reste une forte variabilité d'année en
@@ -438,7 +469,7 @@ moyen est suffisamment représentatif pour permettre d'extraire la composante
 tendancielle lisse des données.
 
 
-## Analyse des variations tendancielles
+### Analyse des variations tendancielles
 
 ```julia; echo=false; output="hidden"
 @info "Fitting underlying trend"
@@ -446,11 +477,14 @@ tendancielle lisse des données.
 
 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.
+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.smooth = data.val .- data.phi_mean;
+data.theta = data.val .- data.phi_mean;
 ```
 
 Même s'il reste des oscillations locales, nous constatons tout de même que la
@@ -459,23 +493,29 @@ de monotonie.
 
 ```julia; echo=false
 plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
-plot!(data.date, data.val,    label="measurements")
-plot!(data.date, data.smooth, label="composante tendancielle")
+plot!(data.date, data.val,   label="mesures")
+plot!(data.date, data.theta, label="theta")
 ```
 
 Nous allons maintenant tenter de caractériser la tendance sous-jacente. Au vu de
-la courbe (concave), nous proposons un modèle quadratique de la forme :
+la courbe (convexe), nous proposons un modèle quadratique de la forme :
 $$
-C(t) = \alpha + \beta t + \gamma t^2
+\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. L'estimation est réalisée sur un
+sous-échantillonnage du jeu de données, comprenant un point (moyen) par an.
+
 ```julia; wrap=false
-model = lm(@formula(smooth ~ date_num + date_num^2),
-            by(data, :year, smooth = :smooth=>mean, date_num = :date_num=>mean))
+sample = by(data, :year, theta = :theta=>mean, date_num = :date_num=>mean)
+model  = GLM.lm(@formula(theta ~ date_num + date_num^2), sample)
 ```
 
-Ce genre de modèle semble correspondre assez bien aux données, avec un assez bon
-niveau de confiance sur l'estimation des paramètres $\alpha$ et $\beta$.
+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
 α,  β,  γ  = coef(model);
@@ -490,51 +530,108 @@ qui se traduit en une augmentation annuelle comprise entre
 `j @printf "%.2f" 365*β1` et `j @printf "%.2f" 365*β2` ppm/an
 avec un niveau de confiance de 95%.
 
-On voit toutefois que l'incertitude sur $\beta$ en particulier est de nature à
-engendrer une perte de prédictibilité du modèle en temps long.
+On voit toutefois que l'incertitude sur les paramètres (en particulier $\beta$)
+est de nature à engendrer une perte de prédictibilité du modèle en temps long.
 
 ```julia; echo=false
-data.smooth_pred = predict(model, data)
+data.theta_pred = predict(model, data)
 
 plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
-plot!(data_raw.date, data_raw.val, label="measurements")
-plot!(data.date, data.smooth_pred, label="smooth model", linewidth=3)
-plot!(data.date, α .+ β1 * data.date_num .+ γ*data.date_num.^2, linecolor="red", label="smooth model (CI 95%)")
+plot!(data_raw.date, data_raw.val, 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%)")
 plot!(data.date, α .+ β2 * data.date_num .+ γ*data.date_num.^2, linecolor="red", label=nothing)
 ```
 
-# Reconstruction du signal complet et prédiction
+### 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
-```julia
-date_num = date2num(Date(lastyear)):10:date2num(today()+Year(5))
+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 calculer $\theta(t)$. On y
+ajoute la fonction de forme annuelle $\phi(t)$ par périodicité.
 
-prediction.smooth = predict(model, prediction)
+```julia; results="raw"
+prediction.theta = predict(model, prediction)
 prediction = join(prediction, avg, on=:day)
 
-prediction.val = prediction.smooth .+ prediction.phi_mean
+prediction.val = 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
 idx = data_raw.date .> Date(lastyear)
 plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
-plot!(data_raw.date[idx], data_raw.val[idx], label="measurements")
-plot!(prediction.date, prediction.smooth,    label="smooth model")
-plot!(prediction.date, prediction.val,       label="predicted value")
+plot!(data_raw.date[idx], data_raw.val[idx], label="mesures")
+plot!(prediction.date, prediction.theta,     label="tendance")
+plot!(prediction.date, prediction.val,       label="prédiction")
 ```
 
-## Pour aller plus loin
+### Pour aller plus loin
+
+```julia; echo=false
+@info "Studying alternate model"
+```
+
+En utilisant des modèles plus complexes, il est possible de mieux caractériser
+la tendance en temps longs :
 
 ```julia; wrap=false
-model = glm(@formula(smooth ~ date_num),
-            by(data, :year, smooth = :smooth=>mean, date_num = :date_num=>mean),
-            InverseGaussian(), InverseSquareLink())
+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
+```julia; echo=false
+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.val = prediction2.theta .+ prediction2.phi_mean
+
+idx = data_raw.date .> Date(lastyear)
+plot(xlabel="Date", ylabel="CO2 [ppm]", legend=:bottomright)
+plot!(data_raw.date[idx], data_raw.val[idx], label="mesures")
+plot!(prediction2.date, prediction2.theta,     label="tendance")
+plot!(prediction2.date, prediction2.val,       label="prédiction")
+```
+
+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).
+
+```julia; echo=false; results="hidden"
+# close plots
 closeall()
 ```
 
diff --git a/module3/exo3/utils.jl b/module3/exo3/utils.jl
new file mode 100644
index 0000000..46b6860
--- /dev/null
+++ b/module3/exo3/utils.jl
@@ -0,0 +1,61 @@
+# Affichage des DataFrames dans la sortie HTML
+function info(df::DataFrame)
+    pretty(x) = x
+    pretty(x::Float64) = @sprintf "%.6f" x
+
+    pretty(::Type{Union{T,Missing}}) where {T} = "$T?"
+
+    function printrow(i)
+        print("")
+        print("$i")
+        for j in 1:ncol(df)
+            print("$(pretty(df[i,j]))")
+        end
+        print("")
+    end
+
+    function printrows()
+        i = 1
+        while i<=nrow(df) && i<=3
+            printrow(i)
+            i += 1
+        end
+
+        i>nrow(df) && return
+
+        if i < nrow(df)-2
+            print("...")
+        end
+
+        i = max(i, nrow(df)-2)
+        while i<=nrow(df)
+            printrow(i)
+            i += 1
+        end
+    end
+
+    println("$(nrow(df)) rows × $(ncol(df)) columns")
+    print("")
+    print("")
+    print("")
+    for col in names(df)
+        print("")
+    end
+    print("")
+    print("")
+    for col in names(df)
+        print("")
+    end
+    print("")
+    print("")
+    print("")
+    printrows()
+    print("")
+    print("
$col
$(pretty(eltype(df[!,col])))
\n")
+end
+
+
+# Permet de construire des durées avec une syntaxe comme :
+# 7Days
+const Years = Year(1)
+const Days  = Day(1)
-- 
2.18.1