Importation des packages

library(tidyr)

Import des données et création de nouvelles variables

# Import des données 
Immatriculations_voitures <- read.csv("h:/Mes documents/Téléchargements/Immatriculations_voitures.csv", sep=";")

# Transformation des données 
Immatriculations_voitures$Date=as.Date(paste(Immatriculations_voitures$Date,"-01", sep=""))  

#-------------- Mois-----------------------------
format(x=Immatriculations_voitures$Date, format="%B") 
##   [1] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##   [7] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [13] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [19] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [25] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [31] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [37] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [43] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [49] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [55] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [61] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [67] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [73] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [79] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [85] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
##  [91] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
##  [97] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [103] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [109] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [115] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [121] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [127] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [133] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [139] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [145] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [151] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [157] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [163] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [169] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [175] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [181] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [187] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [193] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [199] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [205] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [211] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [217] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [223] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [229] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [235] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [241] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [247] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [253] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [259] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [265] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [271] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [277] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [283] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [289] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [295] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [301] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [307] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [313] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [319] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre" 
## [325] "janvier"   "février"   "mars"      "avril"     "mai"       "juin"     
## [331] "juillet"   "août"      "septembre" "octobre"   "novembre"  "décembre"
Immatriculations_voitures$Mois<-format(x=Immatriculations_voitures$Date, format="%B")
Immatriculations_voitures$Mois=as.factor(Immatriculations_voitures$Mois)

#--------------Années ---------------------------------
format(x=Immatriculations_voitures$Date, format="%Y") 
##   [1] "1994" "1994" "1994" "1994" "1994" "1994" "1994" "1994" "1994" "1994"
##  [11] "1994" "1994" "1995" "1995" "1995" "1995" "1995" "1995" "1995" "1995"
##  [21] "1995" "1995" "1995" "1995" "1996" "1996" "1996" "1996" "1996" "1996"
##  [31] "1996" "1996" "1996" "1996" "1996" "1996" "1997" "1997" "1997" "1997"
##  [41] "1997" "1997" "1997" "1997" "1997" "1997" "1997" "1997" "1998" "1998"
##  [51] "1998" "1998" "1998" "1998" "1998" "1998" "1998" "1998" "1998" "1998"
##  [61] "1999" "1999" "1999" "1999" "1999" "1999" "1999" "1999" "1999" "1999"
##  [71] "1999" "1999" "2000" "2000" "2000" "2000" "2000" "2000" "2000" "2000"
##  [81] "2000" "2000" "2000" "2000" "2001" "2001" "2001" "2001" "2001" "2001"
##  [91] "2001" "2001" "2001" "2001" "2001" "2001" "2002" "2002" "2002" "2002"
## [101] "2002" "2002" "2002" "2002" "2002" "2002" "2002" "2002" "2003" "2003"
## [111] "2003" "2003" "2003" "2003" "2003" "2003" "2003" "2003" "2003" "2003"
## [121] "2004" "2004" "2004" "2004" "2004" "2004" "2004" "2004" "2004" "2004"
## [131] "2004" "2004" "2005" "2005" "2005" "2005" "2005" "2005" "2005" "2005"
## [141] "2005" "2005" "2005" "2005" "2006" "2006" "2006" "2006" "2006" "2006"
## [151] "2006" "2006" "2006" "2006" "2006" "2006" "2007" "2007" "2007" "2007"
## [161] "2007" "2007" "2007" "2007" "2007" "2007" "2007" "2007" "2008" "2008"
## [171] "2008" "2008" "2008" "2008" "2008" "2008" "2008" "2008" "2008" "2008"
## [181] "2009" "2009" "2009" "2009" "2009" "2009" "2009" "2009" "2009" "2009"
## [191] "2009" "2009" "2010" "2010" "2010" "2010" "2010" "2010" "2010" "2010"
## [201] "2010" "2010" "2010" "2010" "2011" "2011" "2011" "2011" "2011" "2011"
## [211] "2011" "2011" "2011" "2011" "2011" "2011" "2012" "2012" "2012" "2012"
## [221] "2012" "2012" "2012" "2012" "2012" "2012" "2012" "2012" "2013" "2013"
## [231] "2013" "2013" "2013" "2013" "2013" "2013" "2013" "2013" "2013" "2013"
## [241] "2014" "2014" "2014" "2014" "2014" "2014" "2014" "2014" "2014" "2014"
## [251] "2014" "2014" "2015" "2015" "2015" "2015" "2015" "2015" "2015" "2015"
## [261] "2015" "2015" "2015" "2015" "2016" "2016" "2016" "2016" "2016" "2016"
## [271] "2016" "2016" "2016" "2016" "2016" "2016" "2017" "2017" "2017" "2017"
## [281] "2017" "2017" "2017" "2017" "2017" "2017" "2017" "2017" "2018" "2018"
## [291] "2018" "2018" "2018" "2018" "2018" "2018" "2018" "2018" "2018" "2018"
## [301] "2019" "2019" "2019" "2019" "2019" "2019" "2019" "2019" "2019" "2019"
## [311] "2019" "2019" "2020" "2020" "2020" "2020" "2020" "2020" "2020" "2020"
## [321] "2020" "2020" "2020" "2020" "2021" "2021" "2021" "2021" "2021" "2021"
## [331] "2021" "2021" "2021" "2021" "2021" "2021"
Immatriculations_voitures$Annees<-format(x=Immatriculations_voitures$Date, format="%Y")
Immatriculations_voitures$Annees=as.factor(Immatriculations_voitures$Annees)

# Vérification 
str(Immatriculations_voitures)
## 'data.frame':    336 obs. of  4 variables:
##  $ Date  : Date, format: "1994-01-01" "1994-02-01" ...
##  $ Nombre: int  129170 135855 188301 184343 164836 135987 211104 147451 137127 173721 ...
##  $ Mois  : Factor w/ 12 levels "août","avril",..: 5 4 9 2 8 7 6 1 12 11 ...
##  $ Annees: Factor w/ 28 levels "1994","1995",..: 1 1 1 1 1 1 1 1 1 1 ...
View(Immatriculations_voitures)

Statistiques descriptives

#Stat descriptives
summary(Immatriculations_voitures$Nombre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   20484  148401  168833  167045  186703  292625

Tableau de Buys-Ballot

# Tableau années x mois
Immatriculations_voitures.bb= pivot_wider(data = Immatriculations_voitures, # tableau de données
                     id_cols = Annees, # en ligne
                     names_from = Mois, # en colonnes
                     values_from = Nombre # valeurs
                     )

# affichage
Immatriculations_voitures.bb

Création de la série temporelle pour la période étudiée (1994-2021)

Immat_ts=ts(data = Immatriculations_voitures$Nombre,
              start=c(1994,1),
              frequency=12)

Visualisations & diagnostics graphiques

## visualisation de la série 1994-2021
plot(Immat_ts, col="gray", ylab="Nombre de plaques d'immatriculations crées", 
     xlim=c(1994,2022))
points(Immat_ts, col="blue", pch=20)

#ACF

acf(x = Immat_ts, ylim=c(-1,1),main="ACF de la série")

## Courbes empilées
matplot(t(Immatriculations_voitures.bb[1:27
                                       ,-1]),
        type='b', 
        xlim=c(1,15),
        xlab="Mois", ylab="Nombre de plaques d'immatriculations crées", 
        col=rainbow(27),
        pch=20, lty=1,
        xaxt="n")

axis(side=1,  
     at=1:12,
     labels = format(ISOdate(1994, 1:12, 1), "%b"),
     las=2)

legend("topright",legend =  1994:2021, col=rainbow(27), lty=1)

# moy et ecart-type / année
Immatriculations_voitures.m=apply(X = Immatriculations_voitures.bb[1:27,-1],MARGIN = 2,FUN = mean, na.rm=TRUE)
Immatriculations_voitures.sd=apply(X = Immatriculations_voitures.bb[1:27,-1],MARGIN = 2,FUN = sd, na.rm=TRUE)
# droite des moindres carrés
Immatriculations_voitures.msd.droitemc=lm(Immatriculations_voitures.sd~Immatriculations_voitures.m)

# graph ecart-type en fonction de moy
plot(x=Immatriculations_voitures.m, y=Immatriculations_voitures.sd, 
     pch=20, col="darkblue",  xlab="moyennes par année", ylab="ecarts-types par année")
abline(Immatriculations_voitures.msd.droitemc$coefficients,col="orange",lty=2) # ajout de la droite des moindres carrés
text(x=180000,y=20000,
     labels=paste("coef pente = ", round(Immatriculations_voitures.msd.droitemc$coefficients[2],4)), col="orange")

# Etude de la série approche paramétrique ## Estimation de la tendance

T=length(Immat_ts)
time =1:T
time2 = time^2

# ajustement du modèle polynomial de degré 2: b_0 + b_1*t + b_2*t^2
ajust_tend = lm(formula = Immat_ts~time+time2) # modèle polynomiale degré 2
ajust_tend$coefficients
##   (Intercept)          time         time2 
##  1.607778e+05  1.887165e+02 -6.754424e-01
plot(Immat_ts, col="gray", ylab="Nombre",
     main="ajustement tendance polynomiale (deg 2) / moindres carrés")
points(Immat_ts, col="darkblue", pch=20)

serie_ajustTend = ts(ajust_tend$fitted.values,
              start = c(1994,1),
              frequency=12)
plot(Immat_ts, col="gray", ylab="Nombre",
     main="ajustement tendance polynomiale (deg 2) / moindres carrés")
points(Immatriculations_voitures, col="darkblue", pch=20)
lines(serie_ajustTend, col ="blue", lwd=2)

# série sans tendance
serie_sansTend = Immat_ts - ajust_tend$fitted.values
# Ajustement saisonnalité
S = 12
sint = sin(2*pi*time/S)
cost = cos(2*pi*time/S)

# ajustement du modèle b_0 + b_1*sint + b_2*cost
ajust_sais = lm(serie_sansTend ~ sint + cost)
ajust_sais$coefficients
##   (Intercept)          sint          cost 
## -4.069274e-12  1.042173e+04 -3.594300e+03
## visualisation de la série sans tendance
plot(serie_sansTend, col="gray", ylab="Nombre -- centré", main="ajustement saisonnalité / moindres carrés")
points(serie_sansTend, col="darkblue", pch=20)

## ajout de la tendance estimée
serie_ajustSais = ts(ajust_sais$fitted.values,
                      start = c(1994,1),
                      frequency=12)
lines(serie_ajustSais,col ="blue", lwd=2)

## visualisation de la série
plot(Immat_ts, col="gray", ylab="Nombre ", main="ajustement tendance + saisonnalité / moindres carrés")
points(Immat_ts, col="darkblue", pch=20)

## ajout de la tendance estimée
serie_estim =  serie_ajustTend + serie_ajustSais
lines(serie_estim, col ="blue", lwd=2)

N=48
time_prev=(T+1):(T+N)
serie_2021=window(Immat_ts,
                  start=c(2021,1))
# nouvelles données: les noms des colonnes de ce tableau doivent être *les memes*
# que ceux des modèles
new = data.frame(time=time_prev,
                 time2=time_prev^2,
                 sint=sin(2*pi*time_prev/S),
                 cost=cos(2*pi*time_prev/S)
                 )

# prévision de la tendance
pred_tend= predict(ajust_tend, newdata = new )

# prévision de la saisonnalité
pred_sais= predict(ajust_sais, newdata = new )

# prevision de la série pour les 12 mois de 2021
pred = pred_tend + pred_sais
pred
##        1        2        3        4        5        6        7        8 
## 149764.0 154627.0 157551.9 157682.9 154912.6 149910.7 143944.3 138538.8 
##        9       10       11       12       13       14       15       16 
## 135068.7 134389.8 136609.4 141058.0 146468.4 151315.2 154223.8 154338.6 
##       17       18       19       20       21       22       23       24 
## 151552.1 146534.0 140551.4 135129.6 131643.4 130948.3 133151.7 137584.1 
##       25       26       27       28       29       30       31       32 
## 142978.2 147808.8 150701.2 150799.8 147997.1 142962.8 136964.0 131526.0 
##       33       34       35       36       37       38       39       40 
## 128023.5 127312.2 129499.4 133915.6 139293.5 144107.9 146984.1 147066.5 
##       41       42       43       44       45       46       47       48 
## 144247.6 139197.0 133182.0 127727.8 124209.1 123481.6 125652.6 130052.6
## visualisation de la série
plot(Immat_ts, col="black", ylab="Nombre ",
     xlim=c(1994,2021+4)
     , main="ajustement tendance + saisonnalité / moindres carrés")
points(Immat_ts, col="black", pch=20)

# ajout des prévisions
serie_pred2022 = ts(pred, start=2021, frequency=12)
lines(serie_pred2022, col="red", type='o', pch="*")

# ajout des vraies observations
lines(serie_2021, col="red", type='l', pch=1)
abline(v=2021,lty=3,col="darkgray")

acf(x = Immat_ts, ylim=c(-1,1),main="ACF de la série")

# Suppression de la tendance: 
y=diff(Immat_ts,lag=1)
acf(y,lag.max=40, ylim=c(-1,1), main="")

# Suppression de la saisonnalité: différenciation ordre 12
z=diff(y,lag=12)
acf(z, lag.max=40, ylim=c(-1,1), main="")

# MM ordre 12 // pour retirer la saisonnalité
p=12
MM= stats::filter(x = Immat_ts,
       filter=c(1/(2*p), rep(1/p, p-1), 1/(2*p)))
## visualisation de la série
plot(Immat_ts, col="gray", ylab="Nombre", main="ajustement tendance / moyennes mobiles")
points(Immat_ts, col="darkblue", pch=20)

## ajout de la tendance estimée
lines(MM, col="blue", lwd=2)

##### Avec les moyennes saisonnières
moy_sais=aggregate(serie_sansTend~cycle(serie_sansTend),
                   FUN=mean)

nb_an_obs=27
sais_val=rep(moy_sais$serie_sansTend,nb_an_obs)
sais_serie=ts(sais_val,start=1994, freq=12)

## visualisation de la série sans tendance
plot(serie_sansTend, col="gray", ylab="Nombre -- centré", main="ajustement saisonnalité / moyennes saisonnières")
points(serie_sansTend, col="darkblue", pch=20)

## ajout de la tendance estimée
lines(sais_serie,
      col ="blue", lwd=2)

### Holt-Winters: tendance + saisonnalité
serie_holt=HoltWinters(x = Immat_ts)

pred_holt=predict(serie_holt,
                      n.ahead=N, # 12 mois de 2021
                      prediction.interval = TRUE)
pred_holt
##                fit      upr         lwr
## Jan 2022 104674.39 151726.7  57622.0837
## Feb 2022 121709.21 169552.3  73866.0786
## Mar 2022 139631.00 188252.1  91009.8990
## Apr 2022 102837.81 152224.6  53451.0019
## May 2022 123571.01 173711.8  73430.1783
## Jun 2022 192161.03 243044.7 141277.3509
## Jul 2022 122838.02 174453.9  71222.1736
## Aug 2022  79969.93 132307.7  27632.1639
## Sep 2022 127887.05 180936.9  74837.1945
## Oct 2022 128566.96 182319.5  74814.4446
## Nov 2022 116846.14 171292.3  62400.0253
## Dec 2022 152441.06 207572.0  97310.0776
## Jan 2023 101772.33 161931.1  41613.5182
## Feb 2023 118807.15 179586.5  58027.8094
## Mar 2023 136728.94 198122.5  75335.3366
## Apr 2023  99935.76 161937.5  37933.9781
## May 2023 120668.95 183273.0  58064.9090
## Jun 2023 189258.98 252459.5 126058.4028
## Jul 2023 119935.96 183727.5  56144.4349
## Aug 2023  77067.87 141444.9  12690.8177
## Sep 2023 124984.99 189942.3  60027.6953
## Oct 2023 125664.91 191197.3  60132.4972
## Nov 2023 113944.08 180046.6  47841.5655
## Dec 2023 149539.00 216206.8  82871.2553
## Jan 2024  98870.27 169752.2  27988.3567
## Feb 2024 115905.09 187314.4  44495.7687
## Mar 2024 133826.88 205759.8  61894.0120
## Apr 2024  97033.70 169486.3  24581.0673
## May 2024 117766.90 190735.6  44798.2063
## Jun 2024 186356.92 259838.0 112875.7933
## Jul 2024 117033.90 191023.9  43043.8896
## Aug 2024  74165.81 148661.2   -329.6118
## Sep 2024 122082.94 197080.4  47085.5100
## Oct 2024 122762.85 198258.9  47266.7564
## Nov 2024 111042.02 187033.5  35050.5380
## Dec 2024 146636.95 223120.6  70153.2747
## Jan 2025  95968.21 176151.8  15784.6174
## Feb 2025 113003.04 193653.2  32352.8347
## Mar 2025 130924.82 212038.9  49810.7005
## Apr 2025  94131.64 175707.0  12556.2336
## May 2025 114864.84 196898.9  32830.7420
## Jun 2025 183454.86 265945.1 100964.6247
## Jul 2025 114131.84 197075.7  31187.9763
## Aug 2025  71263.75 154658.8 -12131.2785
## Sep 2025 119180.88 203024.6  35337.1117
## Oct 2025 119860.79 204150.9  35570.6773
## Nov 2025 108139.97 192874.1  23405.8572
## Dec 2025 143734.89 228910.7  58559.0978
## visualisation de la série
plot(Immat_ts, col="gray", ylab="Nombre ",
     xlim=c(1994,2021+4)
     , main="ajustement tendance + saisonnalité et prévisions/ Holt-Winters")
points(Immat_ts, col="darkblue", pch=20)


# ajout de l'estimaiton de la série
lines(serie_holt$fitted[,1], col="blue", type='l', lwd=2)

# ajout des prévisions
lines(pred_holt[,1], col="red", type="o",pch="*")
lines(pred_holt[,2], lty=2, col="orange")
lines(pred_holt[,3], lty=2, col="orange")
abline(v=2022,lty=3,col="darkgray")

# ajout des vraies observations
lines(serie_2021, col="gray", type='o', pch="15")