Analyse Approfondie du Modèle SARIMA
Cheikh Youssou Bamar Gueye
2025-03-23
Introduction
Ce document présente une analyse approfondie de la modélisation SARIMA d’une série temporelle simulée.
Nous explorerons les propriétés de la série, testerons différents modèles, évaluerons leur pertinence et anal-
yserons les résidus.
Simulation des données
Nous générons une série temporelle simulée avec un modèle SARIMA(1,1,0)(1,1,1)[12].
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## [Link] zoo
library(ggplot2)
library(tseries)
library(TSA)
## Registered S3 methods overwritten by 'TSA':
## method from
## [Link] forecast
## [Link] forecast
##
## Attachement du package : 'TSA'
## Les objets suivants sont masqués depuis 'package:stats':
##
## acf, arima
## L'objet suivant est masqué depuis 'package:utils':
##
## tar
library(lmtest)
## Le chargement a nécessité le package : zoo
##
## Attachement du package : 'zoo'
## Les objets suivants sont masqués depuis 'package:base':
##
## [Link], [Link]
1
library(ggfortify)
## Registered S3 methods overwritten by 'ggfortify':
## method from
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
## [Link] forecast
library(seasonal)
library(stats)
library(ggforce)
library(caschrono)
n <- 2000
[Link](123)
# Simulation d'une série SARIMA(1,1,0)(1,1,1)[12]
y <- [Link](n = 2000, list(ar = 0.7, order = c(1, 1, 0)),
seasonal = list(order = c(1, 1, 1), period = 12,
ar = 0.5, ma = 0.4),
innov = rnorm(2000, mean = 0, sd = 1),
[Link] = rnorm(2000, mean = 0, sd = 1))
ts_data <- ts(y, frequency = 12)
length(ts_data)
## [1] 2001
# Affichage de la série
#autoplot(ts_data) + ggtitle("Série temporelle simulée")
Analyse exploratoire
Visualisation de la série
Nous Observons la tendance, la saisonnalité et les fluctuations.
autoplot(ts_data) + ggtitle("Série temporelle simulée")
2
Série temporelle simulée
−50
−100
−150
0 50 100 150
1ère figure : Série temporelle simulée
Cette figure montre une série temporelle générée selon un modèle SARIMA(1,1,0)(1,1,1)[12]. Voici quelques
observations :
• La série semble présenter une tendance globale qui évolue au fil du temps, avec des fluctuations signi-
ficatives.
• On observe des variations irrégulières qui pourraient correspondre à une structure saisonnière ou à des
chocs aléatoires.
• La présence d’une dérive suggère que la série pourrait être non stationnaire, ce qui est cohérent avec
la présence d’une différenciation (d=1) dans le modèle.
Résumé statistique
Nous identifions la Moyenne, la variance, et autres statistiques descriptives.
summary(ts_data)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -174.74 -122.36 -106.85 -102.07 -87.73 14.67
Tendance et Saisonnalite : Décomposition de la série temporelle
Nous utilisons la décomposition STL pour observer les composantes tendance, saisonnière et résiduelle.
decomposed <- stl(ts_data, [Link] = "periodic")
autoplot(decomposed) + ggtitle("Décomposition STL de la série")
3
Décomposition STL de la série
Data
0
−50
−100
−150
remainder
4
0
−4
−8
seasonal
0.1
0.0
−0.1
trend
0
−50
−100
−150
0 50 100 150
2ème figure : Décomposition STL de la série
La décomposition STL permet d’analyser les différentes composantes de la série :
• Données originales ("Data") : Cette première sous-figure affiche la série initiale, qui reprend les
variations observées dans la première figure.
• Composante résiduelle ("Remainder") : Il s’agit de la partie non expliquée par la tendance et la
saisonnalité.
– On observe une série de fluctuations irrégulières, avec des variations d’amplitude relativement
constantes.
• Composante saisonnière ("Seasonal") : Cette courbe met en évidence un motif périodique récur-
rent.
– Le modèle SARIMA(1,1,0)(1,1,1)[12] suppose une périodicité de 12 unités de temps, ce qui est
confirmé ici avec des oscillations régulières.
• Composante tendance ("Trend") : Elle représente l’évolution globale de la série temporelle.
– On observe une tendance descendante initiale suivie d’une stabilisation avec de légères fluctuations,
ce qui indique une dynamique sous-jacente qui évolue dans le temps.
Conclusion
L’analyse de la série temporelle montre qu’elle est influencée par :
• Une forte composante saisonnière (mise en évidence par la décomposition STL).
• Une tendance qui évolue avec le temps.
• Des variations résiduelles qui semblent aléatoires et relativement stationnaires.
4
Test de stationnarité
Nous Confirmons que la série est non stationnaire (tests ADF et KPSS).
[Link](ts_data)
##
## Augmented Dickey-Fuller Test
##
## data: ts_data
## Dickey-Fuller = -3.3135, Lag order = 12, p-value = 0.06833
## alternative hypothesis: stationary
Si on utilise un seuil de significativité de 5 (α = 0.05), la p-value est supérieure à 0.05. Par conséquent, on
ne rejette pas l’hypothèse nulle. Cela signifie que la série n’est pas stationnaire selon le test ADF.
Préparation des données pour la modélisation
Différenciation saisonnière
Nous appliquons une différenciation à l’ordre de la saisonnalité qui est 12 ici.
ts_data_saison = diff(ts_data,lag = 12,differences = 1)
head(ts_data_saison)
## Jan Feb Mar Apr May Jun
## 2 -5.724970 -4.241111 -2.687719 -1.336731 -0.487612 -1.239974
length(ts_data_saison)
## [1] 1989
# Ajuster les marges
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
# Tracer l'ACF avec des couleurs personnalisées
acf(ts_data_saison, main = "ACF de la série différenciée", col = "blue", lwd = 2)
# Tracer la PACF avec des couleurs personnalisées
pacf(ts_data_saison, main = "PACF de la série différenciée", col = "red", lwd = 2)
5
PACF de la série différenciée
ACF de la série différenciée
1.0
1.0
0.8
0.5
0.6
Partial ACF
ACF
0.4
0.0
0.2
−0.5
0.0
0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5
Lag Lag
Une décroissance lente de l’ACF suggère que la série présente une tendance persistante ou une non-
stationnarité. La différenciation saisonnière (D = 1) a déjà été appliquée pour éliminer la tendance. L’ACF
décroît toujours lentement après différenciation, cela pourrait indiquer qu’une deuxième différenciation
non saisonnière est nécessaire (d = 1).Dans le pacf, l’absence de pic significatif à lag 24 suggère que la
saisonnalité n’est pas forte au deuxième cycle, donc une deuxième différentiation saisonnière n’est pas
nécessaire\
• Un pic significatif à lag 12 dans l’ACF confirme la présence d’une saisonnalité annuelle (pour des
données mensuelles).\
• Un pic à lag 12 dans la PACF confirme la présence d’une composante AR saisonnière donc Pmax = 2.\
• Le pic à lag 12 persiste après différenciation, cela pourrait indiquer qu’une composante MA saisonnière
est nécessaire dans le modèle dans l’acf donc Qmax = 1.
Différenciation non saisonnière
Nous appliquons une différenciation standard (d=1) pour éliminer la tendance.
ts_data_dif_saison = diff(ts_data_saison,lag = 1,differences = 1)
head(ts_data_dif_saison)
## Feb Mar Apr May Jun Jul
## 2 1.48385863 1.55339222 1.35098864 0.84911862 -0.75236170 0.03973836
length(ts_data_dif_saison)
## [1] 1988
# Ajuster les marges
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
6
# Tracer l'ACF avec des couleurs personnalisées
acf(ts_data_dif_saison, main = "ACF de la série différenciée", col = "blue", lwd = 2)
# Tracer la PACF avec des couleurs personnalisées
pacf(ts_data_dif_saison, main = "PACF de la série différenciée", col = "red", lwd = 2)
PACF de la série différenciée
ACF de la série différenciée
0.6
0.6
0.4
0.4
Partial ACF
0.2
ACF
0.2
0.0
0.0
−0.2
−0.2
−0.4
0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5
Lag Lag
La décroissance rapide de l’ACF suggère que la série ne présente plus une tendance persistante ou une
stationnarité.\
• Pic significatif à lag 1 : Un pic à lag 1 dans l’acf indique une composante MA(1) donc qmax = 3.\
• Pic significatif à lag 1 : Un pic à lag 1 dans la PACF suggère une composante AR(1) donc pmax = 1.\
• Autres pics : Les pics aux autres lags ne sont pas significatifs, ce qui suggère qu’une composante MA
ou AR d’ordre supérieur n’est pas nécessaire.
Tests de stationnarité :
Nous Confirmons que la série différenciée est stationnaire (tests ADF et KPSS).
[Link](ts_data_dif_saison)
## Warning in [Link](ts_data_dif_saison): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts_data_dif_saison
## Dickey-Fuller = -13.609, Lag order = 12, p-value = 0.01
7
## alternative hypothesis: stationary
Si on utilise un seuil de significativité de 5 (α = 0.05), la p-value est inférieure à 0.05. Par conséquent, on
rejette l’hypothèse nulle. Cela signifie que la série est stationnaire selon le test ADF.
Modélisation SARIMA
Sélection du modèle
Nous Testons plusieurs combinaisons de paramètres et comparer les modèles en utilisant l’AIC et le BIC. Nous
Sélectionnons le modèle avec les valeurs d’AIC et de BIC les plus faibles, tout en évitant le sur-ajustement.
# Ajustez plusieurs modèles ARIMA et SARIMA en utilisant différentes combinaisons d'ordres identifiés
d = 1
D = 1
S = 12
coef <- [Link](pc = c(0:1), qc = c(0:3), Pc = c(0:2), Qc = c(0:1))
coefs <- cbind(coef[,1], d, coef[,2], coef[,3], D, coef[,4])
f_sar <- function(x, coefs, S) {
Arima(x, order = coefs[1:3], seasonal = list(order = coefs[4:6], period = S), [Link] = FALSE, me
}
[Link] <- apply(coefs, 1, f_sar, x = ts_data, S = 12)
# Comparez les critères d'information tels que AIC et BIC pour sélectionner le meilleur modèle
aic <- lapply([Link],function(x) x$aic)%>%unlist
bic <- lapply([Link],function(x) -2*x$loglik + x$nobs*length(x$coef))%>%unlist
par(mfrow = c(1,1))
ac <- order(aic)
plot(aic[ac],type = 'b',pch = 20,axes = F,xlab = '')
axis(1,c(1:length(aic)),paste(coefs[ac,1],coefs[ac,3],coefs[ac,4],coefs[ac,6]),las = 2)
axis(2)
8000
aic[ac]
7000
6000
1001
1201
1301
1111
1311
1221
0301
0321
0211
0101
0121
1220
1320
0220
1110
1310
0120
0110
1100
1300
0021
0300
0100
0010
8
Ajustement du modèle
Nous ajustons le modèle SARIMA avec les paramètres choisis. En fin, nous affichons les résultats (coefficients,
erreurs standards, AIC, BIC).
coefs[[Link](aic),]
## d D
## 1 1 0 0 1 1
coefs[[Link](bic),]
## d D
## 0 1 0 0 1 0
summary([Link][[[Link](aic)]])
## Series: x
## ARIMA(1,1,0)(0,1,1)[12]
##
## Coefficients:
## ar1 sma1
## 0.7070 -1.0000
## s.e. 0.0158 0.0093
##
## sigma^2 = 0.9737: log likelihood = -2824.37
## AIC=5654.74 AICc=5654.75 BIC=5671.52
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02950373 0.9830498 0.7818813 -0.7710447 2.403344 0.09226009
## ACF1
## Training set -0.007137163
[Link][[[Link](aic)]]$coef
## ar1 sma1
## 0.7070357 -0.9999997
[Link][[order(aic)[2]]]$coef
## ar1 ma1 sma1
## 0.71404962 -0.01403087 -0.99999943
Validation du modèle
modfin = [Link][[[Link](aic)]]
t_stat(modfin)
## ar1 sma1
## [Link] 44.61292 -107.9981
## [Link] 0.00000 0.0000
Le coefficient AR(1) est significatif (t-stat = 44.61292, p-value = 0.0000), ce qui indique une forte dépendance
entre les observations successives.\
Le coefficient SMA(1) est également significatif (t-stat = -107.9981, p-value = 0.0000), ce qui indique une
forte composante MA saisonnière.
9
Visualisation des résidus
Nous tracons les résidus pour vérifier qu’ils ressemblent à un bruit blanc.
ggtsdiag(modfin)
Standardized Residuals
4
2
0
−2
0 50 100 150
ACF of Residuals
1.00
0.75
ACF
0.50
0.25
0.00
0 1 2
Lag
p values for Ljung−Box statistic
1.00
p value
0.75
0.50
0.25
0.00
2.5 5.0 7.5 10.0
Lag
Le graphique ggtsdiag montre les résidus standardisés, l’ACF des résidus et les p-values du test de
Ljung-Box.
Les résidus sont bien un bruit blanc, les points dans le graphique des résidus standardisés sont centrés autour
de zéro, l’ACF des résidus montre un pic significatif, et les p-values du test de Ljung-Box sont supérieures
à 0.05
Tests de normalité
Nous vérifions que les résidus suivent une distribution normale (test de Shapiro-Wilk, test de Jarque-Bera).
[Link](modfin$residuals)
##
## Shapiro-Wilk normality test
##
## data: modfin$residuals
## W = 0.99945, p-value = 0.8646
La p-value est supérieure à 0.05, ce qui signifie que les résidus suivent une distribution normale. Cela valide
l’hypothèse de normalité des résidus.
Tests d’autocorrélation
Nous vérifions qu’il n’y a pas d’autocorrélation dans les résidus (test de Ljung-Box).
10
[Link].2(modfin$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)
## Retard p-value
## [1,] 6 0.78104
## [2,] 12 0.88481
## [3,] 18 0.73693
## [4,] 24 0.86625
## [5,] 30 0.91015
## [6,] 36 0.54246
Les p-values pour les retards 6, 12, 18, 24, 30 et 36 sont toutes supérieures à 0.05. Cela signifie qu’il n’y a
pas d’autocorrélation significative dans les résidus aux différents retards testés. Les résidus ressemblent à un
bruit blanc.
Previsions
Découpage des données
train <- window(y, end = 1970)
test <- window(y, start = 1971)
Ajustement du modèle sur les données d’apprentissage
On vérifie que le modèle sur la série tronquée est toujours validé.
mod_train <- Arima(train, order = c(1,1,0),
list(order = c(0,1,1), period = 12),
[Link] = FALSE, method = "CSS-ML")
summary(mod_train)
## Series: train
## ARIMA(1,1,0)(0,1,1)[12]
##
## Coefficients:
## ar1 sma1
## 0.7103 -1.0000
## s.e. 0.0159 0.0095
##
## sigma^2 = 0.9652: log likelihood = -2772.09
## AIC=5550.19 AICc=5550.2 BIC=5566.93
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03184657 0.978676 0.7785944 -0.7818147 2.422851 0.7039879
## ACF1
## Training set -0.004668563
t_stat(mod_train)
## ar1 sma1
## [Link] 44.55324 -104.8612
## [Link] 0.00000 0.0000
[Link].2(mod_train$residuals,nlag=c(6,12,18,24,30,36),type="Ljung-Box",decim=5)
## Retard p-value
11
## [1,] 6 0.71790
## [2,] 12 0.86595
## [3,] 18 0.62207
## [4,] 24 0.80184
## [5,] 30 0.85999
## [6,] 36 0.52126
[Link](mod_train$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod_train$residuals
## W = 0.99945, p-value = 0.8705
Prédictions
pred_model_train = forecast(mod_train,h = 30,level = 95)
summary(pred_model_train)
##
## Forecast method: ARIMA(1,1,0)(0,1,1)[12]
##
## Model Information:
## Series: train
## ARIMA(1,1,0)(0,1,1)[12]
##
## Coefficients:
## ar1 sma1
## 0.7103 -1.0000
## s.e. 0.0159 0.0095
##
## sigma^2 = 0.9652: log likelihood = -2772.09
## AIC=5550.19 AICc=5550.2 BIC=5566.93
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03184657 0.978676 0.7785944 -0.7818147 2.422851 0.7039879
## ACF1
## Training set -0.004668563
##
## Forecasts:
## Point Forecast Lo 95 Hi 95
## 1971 -88.86362 -90.79500 -86.93225
## 1972 -87.30353 -91.12998 -83.47709
## 1973 -86.14783 -91.88721 -80.40846
## 1974 -85.36940 -92.96149 -77.77730
## 1975 -84.82743 -94.17987 -75.47500
## 1976 -84.39585 -95.40642 -73.38528
## 1977 -84.04048 -96.60796 -71.47299
## 1978 -83.83957 -97.86871 -69.81044
## 1979 -83.74557 -99.14907 -68.34206
## 1980 -83.67538 -100.37442 -66.97635
## 1981 -83.58510 -101.50897 -65.66122
## 1982 -83.51244 -102.59799 -64.42688
12
## 1983 -83.61904 -103.81372 -63.42437
## 1984 -83.72920 -104.98480 -62.47360
## 1985 -83.75988 -106.03247 -61.48729
## 1986 -83.82415 -107.07374 -60.57455
## 1987 -83.88076 -108.07100 -59.69053
## 1988 -83.87435 -108.97208 -58.77662
## 1989 -83.82099 -109.79596 -57.84601
## 1990 -83.83460 -110.65914 -57.01006
## 1991 -83.89297 -111.54164 -56.24429
## 1992 -83.93102 -112.38043 -55.48160
## 1993 -83.91761 -113.14614 -54.68907
## 1994 -83.89955 -113.88716 -53.91194
## 1995 -84.04495 -114.77561 -53.31430
## 1996 -84.18265 -115.64069 -52.72462
## 1997 -84.23291 -116.40320 -52.06262
## 1998 -84.31107 -117.17909 -51.44306
## 1999 -84.37756 -117.92943 -50.82570
## 2000 -84.37817 -118.60067 -50.15567
pred_train = pred_model_train$mean
pred_l_train = ts(pred_model_train$lower,start = 1971,frequency = 1)
pred_u_train = ts(pred_model_train$upper,start = 1971,frequency = 1)
Comparaison des prédictions aux observations réelles
[Link](test,pred_train,pred_l_train,pred_u_train,
xlab="t",ylab="sarima",col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topleft",legend=c("X","X_prev"),col=c(1,2,3,3),lty=c(1,1),lwd=c(3,3))
legend("topright",legend=c("int95%_inf","int95%_sup"),col=c(3,3),lty=c(2,2),lwd=c(2,2))
X int95%_inf
−60
X_prev int95%_sup
−80
sarima
−100
−120
1970 1975 1980 1985 1990 1995 2000
13
Erreur de prévision
accuracy(pred_model_train)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03184657 0.978676 0.7785944 -0.7818147 2.422851 0.7039879
## ACF1
## Training set -0.004668563
Conclusion
Le modèle SARIMA sélectionné permet de bien capturer la dynamique de la série. Les tests de validation
montrent que les résidus sont globalement conformes aux hypothèses du modèle, bien que des améliorations
puissent être envisagées.
14