library(tidyr)
# 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)
#Stat descriptives
summary(Immatriculations_voitures$Nombre)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 20484 148401 168833 167045 186703 292625
# 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
Immat_ts=ts(data = Immatriculations_voitures$Nombre,
start=c(1994,1),
frequency=12)
## 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")