Exercice 2.04

Fiche 2
Author

Clément Poupelin

Published

Invalid Date

Modified

March 4, 2025

Intervenant.e.s

Rédaction

Relecture

Setup

Show the code 
# Données
library(dplyr)        # manipulation des données

# Esthétique
library(ggplot2)     ## ggplot
Show the code 
ggTimeSerie <- function(ts, main_title = NULL) {
  df_series <- data.frame(Time = as.numeric(time(ts)), TimeSerie = ts)
  colnames(df_series) <- c("Time", "TimeSerie")
  
  if(is.null(main_title)){
    main <- latex2exp::TeX(paste0("Série $( x_t )_{t=0, ...,n}$ avec n = ", length(ts)))
  } else 
    main <- latex2exp::TeX(paste0(main_title))
  
  p <- ggplot(df_series, aes(x = Time, y = TimeSerie)) +
    geom_line(color = "red") + 
    labs(title = main,
    x = "Time",
    y = "Simulated series") +
    theme_minimal() 
  
  if(length(time(ts(ts))) == length(ts)){
    p <- p
  } else
    p <- p +
    scale_x_continuous(
    breaks = seq(floor(min(df_series$Time)), ceiling(max(df_series$Time)), by = 2),  
    labels = function(x) floor(x)  
  )
  
  return(p)
}
Show the code 
ggACF <- function(ts) {
  acf_data <- acf(ts, plot = FALSE)
  df_acf <- data.frame(Lag = acf_data$lag, ACF = acf_data$acf)
  
  pacf_data <- pacf(ts, plot = FALSE)
  df_pacf <- data.frame(Lag = pacf_data$lag, PACF = pacf_data$acf)  
  
  # Intervalle de confiance
  ci <- qnorm((1 + 0.95) / 2) / sqrt(length(ts))
  
  # ACF 
  p_acf <- ggplot(df_acf, aes(x = Lag, y = ACF)) +
    geom_segment(aes(xend = Lag, yend = 0), color = "red") +
    geom_point(color = "red") +
    labs(title = "Autocorrelation Function (ACF)", x = "Lag", y = "ACF") +
    geom_hline(yintercept = c(-ci, ci), color = "blue", linetype = "dashed") +
    theme_minimal()
  
  # PACF 
  p_pacf <- ggplot(df_pacf, aes(x = Lag, y = PACF)) +  
    geom_segment(aes(xend = Lag, yend = 0), color = "red") +
    geom_point(color = "red") +
    labs(title = "Partial Autocorrelation Function (PACF)", x = "Lag", y = "PACF") +
    geom_hline(yintercept = c(-ci, ci), color = "blue", linetype = "dashed") +
    theme_minimal()
  
  return(list(ACF = p_acf, PACF = p_pacf))
}
Show the code 
Sim_serie <- function(n, a, b, c, d, e, a0 = 0){
  t <- 1:n
  u <- a0 + a*t + b*cos(pi/6*t)  + c*cos(pi/3*t) + d*sin(pi/6*t)  + e*sin(pi/3*t) + rnorm(n) 
  return(ts(u))
}
Show the code 
ggTimeSerie_vs_FittedSerie <- function(ts, fit, main_title = NULL) {

  if (is.null(main_title)) {
    main <- latex2exp::TeX(paste0("Série $( x_t )_{t=0, ...,n}$ vs Série estimée "))
  } else
    main <- latex2exp::TeX(paste0(main_title))
  
  df <- data.frame(Time = seq_along(ts),
                   ts = ts,
                   Fitted = fit$fitted.values)
  colnames(df) <- c("Time", "ts", "Fitted")
  
  p <- ggTimeSerie(df$ts) +
    geom_line(aes(x = df$Time, y = df$Fitted),
              color = "blue",
              linetype = "dashed") +
    scale_color_manual(values = c("red", "blue")) +
    labs(title = main, y = "Series") +
    theme(legend.position = "topleft") +
    annotate(
      "text",
      x = 10,
      y = max(df$ts),
      label = "Série temporelle",
      color = "red"
    ) +
    annotate(
      "text",
      x = 10,
      y = max(df$ts) - 0.2,
      label = "Série estimée",
      color = "blue"
    )
  
  return(p)
}
Show the code 
set.seed(140400)

Données

Le fichier champ.asc est disponible sur le web à l’adresse suivante http://www.math.sciences.univ-nantes.fr/~philippe/lecture/champ.asc mais on aussi été télécharger et sockée dans le dossier Data du repertoire git

expédition mensuelle de champagne en milliers de bouteilles.

Show the code 
# url <- "http://www.math.sciences.univ-nantes.fr/~philippe/lecture/champ.asc"
# champ <- read.csv(url, header = FALSE)
champ <- read.csv("../Data/champ.asc", header = FALSE)

champ.ts <- ts(champ)
Note

sur ce type de série on peut s’attendre à une saisonnalité de 12 pour les 12 mois d’une année.

Comme notre série est de taille 105, si on divise par les 12 mois d’une année, cela veut dire que nous sommes sur une série qui représente les ventes de champagne sur au moins 8 ans.

Analyse exploratoire

on commence par look serie et acf

Show the code 
gridExtra::grid.arrange(ggTimeSerie(champ.ts, main_title = "Serie of Champagne sales "),
                        ggACF(champ.ts)$ACF,
                        ncol=2)

Résultats

D’après les graphiques obtenus ainsi que l’analyse faite à l’exercice précédent, on constate clairement une série de type multplicatif (variance qui explose avec le temps) présentant une tendance linéaire et à priori une saisonnalité de période 12.

On se propose donc de poser dans la fonction ts un l’argument frequency=12 ainsi que de définir une année de début (ici, sans information particulière nous allons choisir un début à l’année 1997 pourillustrer les amnipulations possibles avec la fonction ts).

Show the code 
champ.ts <- ts(champ, frequency = 12, start = 1997)
gridExtra::grid.arrange(ggTimeSerie(champ.ts, main_title = "Serie of Champagne sales "),
                        ggACF(champ.ts)$ACF,
                        ncol=2)

Décomposition multiplicative

Regardons un peu ce qui se passe en transfo log

Par une transformation logarithmique, on se ramène à une décomposition additive. Cette décomposition multiplicative est intéressante lorsqu’on observe une variation linéaire des effets saisonniers

Show the code 
log_champ.ts <- log(champ.ts)

gridExtra::grid.arrange(ggTimeSerie(log_champ.ts, main_title = "Log Serie of Champagne sales "),
                        ggACF(log_champ.ts)$ACF,
                        ncol=2)

Résultats

On peut observé que le passage au log à permis de contrer la croissance en \(t\) de la variance.
En effet, si on pose \(Y_t = t\varepsilon_t\), alors \(Var(Y_t) = t^2Var(\varepsilon_t)\).

Or, avec le passage au log, on aura que \(Var(log(Y_t)) = Var(log(t\varepsilon_t)) = Var(log(t)+log(\varepsilon_t)) = Var(log(\varepsilon_t))\)

Simulation

Pour différentes valeurs des paramètres \((\alpha, \beta, \gamma)\), simuler les séries suivantes de longueur 100 où \((\varepsilon_t)\) est une suite de variables aléatoires i.i.d. \(\mathcal{N}_{(0, 1)}\)


On va donc simuler des séries de la forme suivante \[\begin{align} \alpha t + \beta cos(\frac{2πt}{12}) + \gamma cos(\frac{2πt}{6}) + \beta' cos(\frac{2πt}{12}) + \gamma' cos(\frac{2πt}{6}) + \varepsilon_t \end{align}\]

Note

on utilise Vectorize qui va ???

Show the code 
n <- rep(100,8)
alpha <- rep(c(0.01, 0.05),4)
beta <- rep(c(-1,1,0.1,2),2)
gamma <- rep(c(-0.1,1,2,-0.5),2)
d <- gamma
e <- beta
a0 <- rep(c(0, 7.5, 8, 8.5),2)

Sim_serie.vect <- Vectorize(Sim_serie)
Sim_serie.res <- Sim_serie.vect(n, alpha, beta, gamma, d, e, a0)

On peut déjà constater que cette série à été construite dans l’optique de prendre les périodes que l’on peut détecter avec l’ACF autour du Lag 6 et du Lag 12 et qui sont adéquate à l’aspect de “double pics” présent dans notre série. Le terme \(\alpha\) est de son côté, présent pour prendre en compte la présence de la tendance linéaire.

Comparaison

Comparer l’allure des séries simulées avec la série des ventes de champagne et la série \((log(Ct))\).


Show the code 
plot_list <- list()
for (i in seq_along(alpha)) {
  title_text <- latex2exp::TeX(
    paste(
      "$alpha$ = ",
      alpha[i],
      ", $beta$ = ",
      beta[i],
      ", $gamma$ = ",
      gamma[i],
      "$beta'$ = ",
      d[i],
      ", $gamma'$ = ",
      e[i]
    )
  )
  
  plot_list[[i]] <- ggTimeSerie(Sim_serie.res[, i], title_text)
}

gridExtra::grid.arrange(grobs = plot_list, ncol = 3, nrow = 3)

Show the code 
gridExtra::grid.arrange(ggTimeSerie(champ.ts, main_title = "Serie of Champagne sales "),
                        ggTimeSerie(log_champ.ts, main_title = "Log Serie of Champagne sales "),
                        ncol=2)

Résultats

La série \(log(C_t)\) est la plus adaptée à l’halure des séries simulées car, pour les series simulées, on a pas la variance qui augmente comme pour\(C_t\).

Analyse inférentielle

Sur cette série, calculer les estimateurs de \((\alpha, \beta, \gamma)\) par la méthode des moindres carrés. Que peut-on dire de la qualité du modèle. Peut on modéliser la série des résidus par un bruit blanc?

Show the code 
t <- seq(1, length(log_champ.ts))
mod <- lm(log_champ.ts ~ t + cos((pi/6)*t)  + cos((pi/3)*t) + sin((pi/6)*t)  + sin((pi/3)*t))
mod %>% summary()

Call:
lm(formula = log_champ.ts ~ t + cos((pi/6) * t) + cos((pi/3) * 
    t) + sin((pi/6) * t) + sin((pi/3) * t))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.76814 -0.14722  0.00133  0.18823  0.49000 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      8.1295777  0.0511237 159.018  < 2e-16 ***
t                0.0044624  0.0008376   5.328 6.24e-07 ***
cos((pi/6) * t)  0.3423451  0.0360585   9.494 1.39e-15 ***
cos((pi/3) * t)  0.2575768  0.0358374   7.187 1.26e-10 ***
sin((pi/6) * t) -0.0520224  0.0357329  -1.456    0.149    
sin((pi/3) * t) -0.3516003  0.0358537  -9.807 2.89e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2596 on 99 degrees of freedom
Multiple R-squared:  0.7226,    Adjusted R-squared:  0.7085 
F-statistic: 51.57 on 5 and 99 DF,  p-value: < 2.2e-16

On constate une forte significativité de tout nos termes.

Show the code 
ggTimeSerie_vs_FittedSerie(log_champ.ts, mod)

Résultats

bonne superposition !!!

Conclusion

blablabla on a vu efficacité du passage log pour diminuer variance puis comment procéder pour retrouver un peu une série de ce type

Session info

Show the code 
sessioninfo::session_info(pkgs = "attached")
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 24.04.1 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-03-04
 pandoc   3.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package * version date (UTC) lib source
 dplyr   * 1.1.4   2023-11-17 [1] CRAN (R 4.4.2)
 ggplot2 * 3.5.1   2024-04-23 [1] CRAN (R 4.4.2)

 [1] /home/clement/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────
Show the code 
url = "http://www.math.sciences.univ-nantes.fr/~philippe/lecture/champ.asc"
data = read.csv(url, header = FALSE)

Ct = ts(data, frequency = 12)


QUESTION 1 : Tracer la série \((C_t)\) et sa suite des auto-corrélations empiriques



QUESTION 2 : En utilisant les résultats de l’exercice précédent, peut-on détecter la présence d’une fonction périodique ou d’une tendance dans cette série.


D’après les graphiques obtenus ainsi que l’analyse faite à l’exercice précédent, on constate clairement une série de type multplicatif (variance qui explose avec le temps) présentant une tendance linéaire et une saisonnalité (de période environ 12 peut-être).

QUESTION 3 : Tracer la série \((log(Ct))\) et sa suite des auto-corrélations empiriques


Show the code 
lCt = log(Ct)

On peut observé que le passage au log à permis de contrer la croissance en \(t\) de la variance.
En effet, si on pose \(Y_t = t\varepsilon_t\), alors \(Var(Y_t) = t^2Var(\varepsilon_t)\).
Or, avec le passage au log, on aura que \(Var(log(Y_t)) = Var(log(t\varepsilon_t)) = Var(log(t)+log(\varepsilon_t)) = Var(log(\varepsilon_t))\)


QUESTION 4 : Pour différentes valeurs des paramètres \((\alpha, \beta, \gamma)\), simuler les séries suivantes de longueur 100 où \((\varepsilon_t)\) est une suite de variables aléatoires i.i.d. \(\mathcal{N}_{(0, 1)}\)


On va donc simuler des séries de la forme suivante \[\begin{align} \alpha t + \beta cos(\frac{2πt}{12}) + \gamma cos(\frac{2πt}{6}) + \beta' cos(\frac{2πt}{12}) + \gamma' cos(\frac{2πt}{6}) + \varepsilon_t \end{align}\]

Show the code 
ut = function(n, a, b, c, d, e, a0 = 0){
  t = 1:n
  eps_t = rnorm(n)
  
  u = a0 + a*t + b*cos(pi/6*t)  + c*cos(pi/3*t) + d*sin(pi/6*t)  + e*sin(pi/3*t) + eps_t 
  return(u)
}

# Simulations pour différentes valeurs de coefficients
n = rep(100,8)
alpha = rep(c(0.01, 0.05),4)
beta = rep(c(-1,1,0.1,2),2)
gamma = rep(c(-0.1,1,2,-0.5),2)
d = gamma
e = beta
a0 = rep(c(0, 7.5, 8, 8.5),2)


# On stocke les simulations 
ut.vect = Vectorize(ut)
simu.res = ut.vect(n,alpha,beta,gamma,d,e,a0)

On peut déjà constater que cette série à été construite dans l’optique de prendre les périodes que l’on peut détecter avec l’ACF autour du Lag 6 et du Lag 12 et qui sont adéquate à l’aspect de “double pics” présent dans notre série. Le terme \(\alpha\) est de son côté, présent pour prendre en compte la présence de la tendance linéaire.


QUESTION 5 : Comparer l’allure des séries simulées avec la série des ventes de champagne et la série \((log(Ct))\).



QUESTION 6 : Pour laquelle des deux séries \(((Ct))\) \((log(Ct))\), le modèle défini en question 4 vous semble le plus pertinent.


La série \(log(C_t)\) est la plus adaptée à l’halure des séries simulées car, pour les series simulées, on a pas la variance qui augmente comme pour\(C_t\).

QUESTION 7 : Sur cette série, calculer les estimateurs de \((\alpha, \beta, \gamma)\) par la méthode des moindres carrés. Que peut-on dire de la qualité du modèle. Peut on modéliser la série des résidus par un bruit blanc?


Show the code 
#construisons un modèle de regression
t = seq(1,length(Ct))

model = lm(log(Ct) ~ t + cos((pi/6)*t)  + cos((pi/3)*t) + sin((pi/6)*t)  + sin((pi/3)*t))
summary(model)

Call:
lm(formula = log(Ct) ~ t + cos((pi/6) * t) + cos((pi/3) * t) + 
    sin((pi/6) * t) + sin((pi/3) * t))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.76814 -0.14722  0.00133  0.18823  0.49000 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      8.1295777  0.0511237 159.018  < 2e-16 ***
t                0.0044624  0.0008376   5.328 6.24e-07 ***
cos((pi/6) * t)  0.3423451  0.0360585   9.494 1.39e-15 ***
cos((pi/3) * t)  0.2575768  0.0358374   7.187 1.26e-10 ***
sin((pi/6) * t) -0.0520224  0.0357329  -1.456    0.149    
sin((pi/3) * t) -0.3516003  0.0358537  -9.807 2.89e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2596 on 99 degrees of freedom
Multiple R-squared:  0.7226,    Adjusted R-squared:  0.7085 
F-statistic: 51.57 on 5 and 99 DF,  p-value: < 2.2e-16

On constate une forte significativité de tout nos termes.

Show the code 
# Pour un complément d'information, on peut également regarder la PACF de nos séries $C_t$ et $log(C_t)$
par(mfrow=c(1, 2))  
pacf(Ct, main = TeX("PACF for serie $(C_t)$"), col = "cyan2")
pacf(lCt, main = TeX("PACF for serie $log(C_t)$"), col = "cyan4")

Show the code 
par(mfrow=c(1, 1))