---
title: "Agen"
author: "Nous"
format: html
editor: visual
---

## Températures records à Agen

Commençons par importer les données **après avoir lu le format des données etc...**.

```{r}
rm(list=ls())
url <- "https://mribatet.perso.math.cnrs.fr/docs/EVT/TG_SOUID260218.txt"
data <- read.table(url, skip = 20,
                   header = TRUE,
                   sep = ",", 
                   na.strings = "-9999")
data$TG <- data$TG / 10##pour être en degrés Celsius
```

Conseil : on regarde si les 1ères et dernières lignes sont bien importées.

Gérons les dates:

```{r}
data$DATE <- as.Date(as.character(data$DATE), 
                     format = "%Y%m%d") 
```

## Analyse descriptive (à l'arrache)

```{r}
summary(data)
```

Il a tout traité en variables quantitatives. Or seul TG l'est réellement.

Il y a 4 valeurs manquantes sinon c'est OK et toutes les années sont complètes donc on va faire comme si tout était nickel !

## Analyse GEV

Construisons les blocs maxima sur les trimestres (exprès pour que le modèle GEV stationnaire soit pourri)

```{r}
library(lubridate)
library(evd)
trim_max <- aggregate(TG ~  year(DATE)*quarters(DATE),
                      FUN = max, data = data)
colnames(trim_max) <- c("Years", "Quarters", "TG")
```

On peut regarder l'évolution des records par trimestres
Ajustons notre GEV (stationnaire) sur ces données
```{r}
boxplot(TG ~ Quarters, data = trim_max)
```
On voit, comme attendu, un effet du trimestre que le modèle purement stationnaire ci dessous ne pourra jamais reproduire.

```{r}
fit_stat <- fgev(trim_max$TG)
fit_stat
```

Adéquation du modèle

```{r}
plot(fit_stat)
```

Le modèle est bien bien nul (cf. qq-plot, pp-plot). Tentons un modèle non stationnaire. Pour être plus précis, je vais faire un modèle où les paramètres d'échelle $\sigma$ et de forme $\xi$ de la GEV sont constants mais où le paramètre de position varie au cours du temps selon la forme suivante: 

$$\mu(t) = \beta_0 + \beta_1 1_{\{t \in \{Q_2, Q_3\}\}} + \beta_2 1_{\{t \in Q_4\}}$$

```{r}
## Construction de la matrice de design (sans la 1ère colone)
design_mat <- data.frame(Q23 = trim_max$Quarters %in% c("Q2", "Q3"),
                         Q4 = trim_max$Quarters == "Q4")
fit_nonstat <- fgev(trim_max$TG, nsloc = design_mat)
fit_nonstat
```
Adéquation du modèle
```{r}
plot(fit_nonstat)
```

Le modèle semble correct.

Prenons un modèle un plus compliqué encore où
$$\mu(t) = \beta_0 + \beta_1 1_{\{t \in \{Q_2, Q_3\}\}} + \beta_2 1_{\{t \in Q_4\}} + \beta_3 \times \frac{t - 1986.5}{100} + \beta_4 \times \left(\frac{t - 1986.5}{100}\right)^2$$

```{r}
design_mat2 <- cbind(trend = (trim_max$Years - 1986.5) / 100, 
                     trend2 = ((trim_max$Years - 1986.5) / 100)^2,
                     design_mat)
fit_nonstat2 <- fgev(trim_max$TG, nsloc = design_mat2)
fit_nonstat2
```

Ce modèle plus compliqué est il adéquat ?
```{r}
plot(fit_nonstat2)
```

Oui ! Il est même étonnement bon.

Y-a-t-il un effet du changement climatique à Agen ? On peut tester cela via un test du rapport de vraisemblance.
```{r}
anova(fit_nonstat2, fit_nonstat)
```

Ici, la p-valeur est inférieur au niveau standard 5% et l'on rejette $H_0$ et l'on garde donc le modèle le plus complexe. Il semblerait bien qu'il y ait un effet du changement climatique sur les températures records à Agen.

On peut même tracer l'évolution du niveau de retour de période de retour 50 ans au cours du temps.

```{r}
ret_levQ23 <- function(year){
  param <- fit_nonstat2$param
  loc <- param["loc"] + param["loctrend"] * (year - 1986.5) / 100 +
    param["loctrend2"] * ((year - 1986.5) / 100 )^2 + param["locQ23"]
  scale <- param["scale"]
  shape <- param["shape"]
  qgev(1 - 1 / (ret_per * 4), loc, scale, shape)## x 4 car 4 trimestres dans une année ;-)
}

ret_per <- 50
plot(ret_levQ23, from = 1949, to = 2050, ylim = c(0, 40),
     xlab = "Années", ylab = "Températures (C)")
tmp <- subset(trim_max, Quarters %in% c("Q2", "Q3"))
points(tmp$Years, tmp$TG)
```
Devoir à la maison pour s'entrainer :
  + Le terme de tendance quadratique est-il nécessaire ?
  + Est ce qu'un modèle de Gumbel non stationnaire serait approprié ici ?


## Rien à voir : Intervalle de confiance

On se fait un intervalle de confiance pour l'espérance. On suppose que $X_1, \ldots, X_n$ sont des copies indépendantes d'une $N(\mu, 1)$. On sait qu'un IC symétrique pour $\mu$ est de la forme
$$\left[\bar{X} - z_{\alpha} \times \frac{1}{\sqrt{n}}, \bar{X} + z_{\alpha} \times \frac{1}{\sqrt{n}}\right],$$
où $z_\alpha$ est un quantile de la loi N(0,1).

Faisons l'expérience numérique une fois
```{r}
mu_star <- -13
nobs <- 50
data <- rnorm(nobs, mu_star, 1)
ic <- mean(data) + c(-1, 1) * qnorm(0.975) / sqrt(nobs)

print(ic)
appartient <- mu_star > ic[1] & mu_star < ic[2]
print(appartient)
```
Notons qu'à chaque expérience, l'IC change !

Refaisons cette expérience un grand nombre de fois et calculons la proportions où $\mu_*$ appartient à l'IC.
```{r}
mu_star <- -13
nobs <- 50
nexp <- 999

appartient <- rep(NA, nexp)
all_ic <- matrix(NA, nexp, 2)
for (i in 1:nexp){
  data <- rnorm(nobs, mu_star, 1)
  ic <- mean(data) + c(-1, 1) * qnorm(0.975) / sqrt(nobs)
  all_ic[i,] <- ic
  appartient[i] <- mu_star > ic[1] & mu_star < ic[2]
}

col <- c("red", "lightgrey")[1 + appartient]
plot(0, xlim = range(all_ic), ylim = c(1, nexp),
     xlab = expression(mu), ylab = "Experience number")
arrows(all_ic[,1], 1:nexp, all_ic[,2], 1:nexp,
       angle = 90, code = 3, col = col, length = 0.05)
abline(v = mu_star, lwd = 2)
print(mean(appartient))
```

