Projet Forecasting
BROSSAY Arnaud & HOPSORE Théo
Introduction
Ce projet est divisé en plusieurs fichiés :
Le fichier que vous consultez actuellement, le raport, faisant une analyse pas à pas de notre code complet.
Un fichier .rmd contenant le code sous forme de notebook R, le raport est donc une version pdf de ce dernier
Un fichier R contenant le code brut
Ce rapport sera divisé en deux sous-partie :
Une partie simulation qui a pour but de construire de toute pièce un modèle ARMA en nous basant sur les connaissances acquises
Une partie pratique qui consiste à analyser et déduire un modèle à partir d’un jeu de données
Partie 1 : Simulation
Nous allons simuler ici un processus \(ARMA(2,2)\) qui peut s’écrire sous la forme suivante : \[X_t-\phi_1X_{t-1}-\phi_2X_{t-2} = \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2}\]
# Charger la bibliothèque de séries chronologiques (time series)
library(stats)
# Définir les paramètres ARMA(2,2)
phi <- c(0.7, -0.3) # Coefficients AR
theta <- c(0.4, -0.2) # Coefficients MA
# Nombre d'observations dans la série temporelle
n <- 100
# Générer la série temporelle (ici, nous utilisons une série ARMA(2,2) pour l'exemple)
set.seed(123)
epsilon <- rnorm(n)
y <- numeric(n)
# Initialiser les premières valeurs (peuvent être aléatoires)
y[1] <- rnorm(1)
y[2] <- rnorm(1)
# Générer les observations ARMA(2,2)
for (i in 3:n) {
y[i] <- phi[1] * y[i-1] + phi[2] * y[i-2] + epsilon[i] + theta[1] * epsilon[i-1] + theta[2] * epsilon[i-2]
}
# Créer un objet de série chronologique
arma_series <- ts(y)
plot(arma_series, type = "l", main = "Simulation d'un processus ARMA(2,2)", ylab = "Valeurs")
acf(arma_series)
pacf(arma_series)
define_future_values_ARMA <- function(data, phi, theta, p, q, h) {
n <- length(data)
future_values <- numeric(h)
for (i in 1:h) {
if (i <= max(p, q)) {
prediction <- sum(phi[1:p] * data[(n-i+1):(n-1)]) + sum(theta[1:q] * rnorm(q))
} else {
prediction <- sum(phi * future_values[(i-p):(i-1)]) + sum(theta * rnorm(q))
}
future_values[i] <- prediction
}
return(future_values)
}future = define_future_values_ARMA(arma_series, phi, theta, 2, 2, 10)# Plot historical data and future predictions
plot(arma_series, type = "l", main = "Historical Data and ARMA(2,2) Predictions", ylab = "Values", xlab = "Time")
lines(length(arma_series):(length(arma_series) + length(future) - 1), future, col = "red")
legend("topright", legend = c("Historical Data", "Predictions"), col = c("black", "red"), lty = 1)
arma_simulated <- arima.sim(model = list(ar = phi, ma = theta), n = n)
# Créer un objet de série chronologique
set.seed(123)
arma_simulated <- ts(arma_simulated)
# Afficher le graphique de la série temporelle
plot(arma_simulated, type = "l", main = "Simulation d'un processus ARMA(2,2)", ylab = "valeurs")
Partie 2 : Forecasting Pratique
Nous commencons tout d’abord par initialiser les librairies essentielles :
timeSeries : donne de nombreux outils de manipulation pour les séries temporelles. Il va notement nous permettre de former un objet série temporelle à l’aide d’un fichier .csv
forecast : ropose des méthodes et des outils pour afficher et analyser des prévisions de séries temporelles.
ggplot2 : permet de visualiser un jeu de données ou des modèles mathématiques.
lmtest : permet de tester la significativité de nos coefficients des modèles ARMA, ARIMA ou SARIMA
tseries : fonction pour tester la stationnarité des séries temporelles
library(timeSeries)
library(forecast)
library(dplyr)
library(ggplot2)
library(tseries)
library(lmtest)Présentation de nos données.
Nous allons nous baser sur le dataset Industrial production of electric and gas utilities in the United States, qui représente l’index de production industrielle en gaz et en électricité sur le territoire américain de manière mensuelle dès 1939.
Le fichier peut être retrouvé sur le site de la Federal Reserve Bank of St.Louis en cliquant ici.
Nous allons faire un premier plot de la série temporelle pour essayer de voir si cette derniere comporterait potentiellement une saisonnalite et/ou une tendance.
data = read.csv('electricity.csv') # lecture du fichier csv
datanames(data)[names(data) == "IPG2211A2N"] <- "Production_Enérgétique"
#on renomme pour plus de lisibilitéPour ce faire, après avoir lu le fichier .csv, nous allons créer un objet timeSeries à partir de ce dernier, nous considérons l’entièreté des observations, donc de 1939 à 2023, observations faites sur une base mensuelle d’ou la frequency à 12.
#passage du fichier csv en série temporelle
data <- ts(
data$Production_Enérgétique,
start = 1939,
end = 2023,
frequency = 12
)
X <- dataplot(X, type = "l", xlab = "Date", ylab = "Production_Enérgétique", main = "Production énergétique entre 1940 et 2022")
Pour faciliter notre analyse et avoir un modèle plus précis, nous n’allons utiliser que les observations à partir de 1985.
X = tail(X, n=457)# on prend les 456 derniers mois
X = head(X,n=456)# pour supprimer le premier mois de décembre, la dernière date que l'on prédira plus tard
plot(X, type = "l", xlab = "Date", ylab = "Production_Enérgétique", main = "Production énergétique entre 1985 et 2022")
On peut supposer à premiere vue que le modèle assimilé à cette série n’est pas stationnaire, en effet, la variance semble augmenter au cours du temps et on observe aussi une tendance à la hausse.
Nous allons ensuite effectuer une décomposition de la série afin de pouvoir mieux la cerner. Ainsi nous allons utiliser la fonction decompose() afin de pouvoir la décomposer en fonction de:
sa tendance
sa saisonalité
son bruit
decomp<-decompose(X)# décomposition en trend,saisonnalité et noise
plot(decomp)
On observe donc que cette série temporelle a une tendance croissante légérement supérieur à une tendance linéaire (trend) et on peut clairement aussi voir la partie saisonnière des données (seasonal)
Nous sommes sur un modèle multipliactif.Celui-ci entraine un accroisement de la variabilité, ce qui explique la transformation logarithmique opérée par la suite. De manière générale, il est préférable de passer au log pour travailler sur la série temporelle.
Y = log(X) #passage au log de notre série car modèle multiplicatif
plot(Y, type = "l", xlab = "Date", ylab = "Production_Enérgétique", main = "ln(X)")
On voit que la série transformée présente toujours une tendance et conserve une saisonnalité (de période 12). On doit donc effectuer une différentiation pour rendre cette série stationnaire.
plot(Acf(Y,lag.max=36,plot=FALSE),lwd=5,type="h",lend="butt",main = "ln(X)")
La sortie ACF présente une décroissance lentre vers 0, ce qui traduit une tendance et donc un problème de non-stationnarité. On effectue une différenciation \((I - B)\) afin de supprimer cette tendance.
#Série différencier une seule fois
y_dif1=diff(Y,lag=1,differences=1)
#lag=1: retard, differences=1: pour indiquer que l’on fait une seule fois
plot(Acf(y_dif1,lag.max=36,plot=FALSE),lwd=5,type="h",lend="butt",main = "(I - B)ln(X)")
Nous pouvons observer des pics positifs importants tous les 6 mois (6,12,18…) ainsi que des pics négatifs important tous les 6 mois en décalage (3,9,15…), nous pensons donc faire face à une saisonalité de 12.
On devrait donc appliquer une différentiation à l’ordre 12 afin de supprimer cette saisonalité.
Différentation \((1-B^{12})Y_t\)
#Série différencier une seconde fois
y_dif_1_12=diff(y_dif1,lag=12,differences=1)
#lag=12: retard
plot(Acf(y_dif_1_12,lag.max=72,plot=FALSE),lwd=5,type="h",lend="butt",main = "(I - B^12)ln(X)")
plot(Pacf(y_dif_1_12,lag.max=72,plot=FALSE),lwd=5,type="h",lend="butt",main = "(I - B^12)ln(X)")
Nous pouvons de nouveau observer des barres se dégager du lot sur les graphes ACF et PACF.
Dans un cas comme dans l’autres, les observations qui resortent se situent sur un lag h multiple de 12 (1, 12, 24…)
En nous basant sur ces observations, nous pouvons imaginer un premier modèle pour représenter notre série : le \(SARIMA(1,1,1)(1,1,1)_{12}\)
En effet, nous avons du effectuer une premiere différentiation de lag 1 pour supprimer la tendance, puis une seconde de lag 12 pour supprimer la saisonalité.
Etude du modèle \(SARIMA(1,1,1)(1,1,1)_{12}\)
model1=arima(Y,order=c(1,1,1),list(order=c(1,1,1),period=12),include.mean=FALSE,method="CSS-ML") #création du premier modèle SAMIRA(1,1,1)(1,1,1)-12
summary(model1) # affiache des coeffs##
## Call:
## arima(x = Y, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12),
## include.mean = FALSE, method = "CSS-ML")
##
## Coefficients:
## ar1 ma1 sar1 sma1
## 0.4669 -0.9164 0.0820 -0.8506
## s.e. 0.0540 0.0262 0.0557 0.0290
##
## sigma^2 estimated as 0.0006392: log likelihood = 992.91, aic = -1975.82
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001653991 0.02492917 0.01949609 -0.03845597 0.4317816 0.2589603
## ACF1
## Training set -0.01237824On obtient différents coefficients ar1, ma1, sar1, sma1, les coefficients avec les “s” devant sont ceux de la partie saisonnière.
Test de significativité des paramètres du model1
test_results <- coeftest(model1) # test des coeffs
test_results##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.466936 0.054042 8.6402 <2e-16 ***
## ma1 -0.916437 0.026223 -34.9483 <2e-16 ***
## sar1 0.082000 0.055662 1.4732 0.1407
## sma1 -0.850645 0.029046 -29.2862 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1On constate que le coefficient sar1 n’est pas significative. En effet sa p-valeur n’est pas inférieur à 5%.
Test de la blancheur des résidus de ce modèle
x1<-model1$residuals #test de la blancheur des résidus du premier modèle
nlag=c(6,12,18,24,30,36)
res = c()
for (i in nlag){
a1=Box.test(x1, lag=i, type = "Ljung-Box")
res=rbind(res,t(rbind(i,a1)[,3]))
}
colnames(res)=c("nlag","p-value x1")
res## nlag p-value x1
## [1,] 6 0.8398823
## [2,] 12 0.6102496
## [3,] 18 0.7003898
## [4,] 24 0.1009543
## [5,] 30 0.07573895
## [6,] 36 0.1643095Nous obtenons ici les p-valeurs pour les ordres 6, 12, 18, ...
• les p-valeurs ne sont pas inférieurs au seuil de 5%, ainsi on ne rejette pas l’hypothèse nulle.
Ce modèle était satisfaisant pour les résidus mais pas pour les paramètres, certains n’etaient pas significatifs.
Ce qu’on doit donc faire :
- On retire le paramètre le moins significatif (le coefficients saisonnier auto-régressif d’odre 1 ) et on ré-étudiie le modèle ainsi obtenu. D’ou l’introduction du second modèle.
Modèle 2 : \(SARIMA(1,1,1)(0,1,1)_{12}\)
model2=arima(Y,order=c(1,1,1),list(order=c(0,1,1),period=12),
#modèle2 SARIMA(1,1,1)(0,1,1)-12
include.mean=FALSE,method="CSS-ML")
summary(model2)##
## Call:
## arima(x = Y, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),
## include.mean = FALSE, method = "CSS-ML")
##
## Coefficients:
## ar1 ma1 sma1
## 0.4612 -0.9142 -0.8309
## s.e. 0.0543 0.0266 0.0279
##
## sigma^2 estimated as 0.0006419: log likelihood = 991.81, aic = -1975.63
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001616266 0.02498259 0.01951662 -0.03770296 0.4322097 0.259233
## ACF1
## Training set -0.01127056test_results2 <- coeftest(model2)
test_results2##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.461216 0.054274 8.4979 < 2.2e-16 ***
## ma1 -0.914196 0.026564 -34.4148 < 2.2e-16 ***
## sma1 -0.830883 0.027869 -29.8142 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1On observe que tous les coefficients sont mainteant significatifs avec un résultat de p-value au test de significativité des coef largement inférieur à 5%, on peut alors maintenant regarder les résidus.
x2<-model2$residuals #test résidu du second modèle
nlag=c(6,12,18,24,30,36)
res = c()
for (i in nlag){
a1=Box.test(x2, lag=i, type = "Ljung-Box")
res=rbind(res,t(rbind(i,a1)[,3]))
}
colnames(res)=c("nlag","p-value x2")
res## nlag p-value x2
## [1,] 6 0.8442835
## [2,] 12 0.4848878
## [3,] 18 0.5610006
## [4,] 24 0.04112886
## [5,] 30 0.02635897
## [6,] 36 0.0593968Malheuresement certaines p-value sont inférieurs au seul de 5%, on rejette donc l’hypothèse nulle et les résidus ne sont pas blancs (c’est à dire qu’ils présentent de l’autocorrélation). Cela suggère que le modèle ne parvient pas à expliquer toute la variabilité de la série temporelle.
On doit donc trouver un autre modèle. Pour cela nous allons repartir de notre premier modèle et nous allons tester plusieurs alternatives.
Modèle 3 : \(SARIMA(2,1,1)(1,1,1)_{12}\)
model3=arima(Y,order=c(2,1,1),list(order=c(1,1,1),period=12),include.mean=FALSE,method="CSS-ML")
#modèle3 SARIMA(2,1,1)(1,1,1)-12
summary(model3)##
## Call:
## arima(x = Y, order = c(2, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12),
## include.mean = FALSE, method = "CSS-ML")
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.4633 0.0192 -0.9206 0.0827 -0.8496
## s.e. 0.0541 0.0525 0.0271 0.0557 0.0292
##
## sigma^2 estimated as 0.0006391: log likelihood = 992.98, aic = -1973.96
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001688407 0.02492663 0.01948926 -0.03922388 0.431655 0.2588695
## ACF1
## Training set -0.005673629test_results3 <- coeftest(model3)
test_results3##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.463305 0.054090 8.5654 <2e-16 ***
## ar2 0.019225 0.052505 0.3662 0.7143
## ma1 -0.920588 0.027068 -34.0108 <2e-16 ***
## sar1 0.082676 0.055695 1.4845 0.1377
## sma1 -0.849574 0.029210 -29.0852 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Les coeffs ar2 et sar1 sont peu significatifs alors on peut essayer autres chose
Modèle 4 : \(SARIMA(1,1,2)(1,1,1)_{12}\)
model4=arima(Y,order=c(1,1,2),list(order=c(1,1,1),period=12),
include.mean=FALSE,method="CSS-ML")
#modèle4 SARIMA(1,1,2)(1,1,1)-12
summary(model4)##
## Call:
## arima(x = Y, order = c(1, 1, 2), seasonal = list(order = c(1, 1, 1), period = 12),
## include.mean = FALSE, method = "CSS-ML")
##
## Coefficients:
## ar1 ma1 ma2 sar1 sma1
## 0.5095 -0.9678 0.0428 0.0829 -0.8493
## s.e. 0.1171 0.1318 0.1089 0.0557 0.0293
##
## sigma^2 estimated as 0.0006391: log likelihood = 992.99, aic = -1973.98
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001693235 0.02492643 0.01948859 -0.0393324 0.4316439 0.2588606
## ACF1
## Training set -0.004800352test_results4 <- coeftest(model4)
test_results4##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.509492 0.117077 4.3518 1.350e-05 ***
## ma1 -0.967796 0.131832 -7.3412 2.118e-13 ***
## ma2 0.042812 0.108909 0.3931 0.6942
## sar1 0.082889 0.055724 1.4875 0.1369
## sma1 -0.849346 0.029292 -28.9960 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1même chose que pour le modèle précédent, non significatif
Modèle 5 : \(SARIMA(1,1,1)(1,1,2)_{12}\)
model5=arima(Y,order=c(1,1,1),list(order=c(1,1,2),period=12),
include.mean=FALSE,method="CSS-ML")
#modèle5 SARIMA(1,1,1)(1,1,2)-12
summary(model5)##
## Call:
## arima(x = Y, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 2), period = 12),
## include.mean = FALSE, method = "CSS-ML")
##
## Coefficients:
## ar1 ma1 sar1 sma1 sma2
## 0.4658 -0.9176 -0.3989 -0.3262 -0.4454
## s.e. 0.0534 0.0254 0.1893 0.1767 0.1433
##
## sigma^2 estimated as 0.0006318: log likelihood = 995.44, aic = -1978.87
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001693182 0.02478505 0.01942303 -0.03931007 0.4301974 0.2579899
## ACF1
## Training set -0.01190318test_results5 <- coeftest(model5)
test_results5##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.465841 0.053352 8.7314 < 2.2e-16 ***
## ma1 -0.917634 0.025408 -36.1158 < 2.2e-16 ***
## sar1 -0.398903 0.189300 -2.1073 0.035095 *
## sma1 -0.326224 0.176746 -1.8457 0.064932 .
## sma2 -0.445371 0.143331 -3.1073 0.001888 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sur ce modèle on a une AIC la plus basse et les coefficients sont tous assez significatif, on peut donc essayer d’étudier les résidus.
Test de blancheur des résidus du modèle 5, \(SARIMA(1,1,1)(1,1,2)_{12}\)
x5<-model5$residuals #test de blacnheur du modèle 5, SARIMA(1,1,1)(1,1,2)-12
nlag=c(6,12,18,24,30,36)
res = c()
for (i in nlag){
a1=Box.test(x5, lag=i, type = "Ljung-Box")
res=rbind(res,t(rbind(i,a1)[,3]))
}
colnames(res)=c("nlag","p-value x2")
res## nlag p-value x2
## [1,] 6 0.8459775
## [2,] 12 0.593727
## [3,] 18 0.7125789
## [4,] 24 0.4003716
## [5,] 30 0.3091912
## [6,] 36 0.4372095Toute les p-value sont supérieur a 5% donc dans ce modèle, le terme d’erreur est bien un bruit blanc.
Vérification de la normalité des résidus
shapiro.test(model5$residuals) # test de la normalité des résidus##
## Shapiro-Wilk normality test
##
## data: model5$residuals
## W = 0.99511, p-value = 0.16La p-valeur est supérieure à 5%, donc on considère ici que le test de normalité n’est pas rejeté. Le résidu est potentiellement gaussien.
Choix définitif d’un modèle en donnant la forme explicite de son équation
Le modèle \(SARIMA(1,1,1)(1,1,2)_{12}\) est le modèle nous donnant la plus basse AIC et dont les résidus valide nos tests (test de blancheur et vérification de la normalité des résidus). Nous allons donc choisir ce modèle.
Sa forme explicite est la suivante :
\((1−ϕ_1B)(1−Φ_1B^{12})(1−B)(1−B^{12})yt=(1+θ_1B)(1+Θ_1B^{12})ϵt\)
Prévision à l’aide du modèle retenu
pred_model5=forecast(model5,h=12,level=95) #prévision avec le modèle 5
pred=exp(pred_model5$mean) #valeur prédite, passage à exp pour inverser le log
pred_l=ts(exp(pred_model5$lower),start=c(2023,1),frequency=12) #estimation haute
pred_u=ts(exp(pred_model5$upper),start=c(2023,1),frequency=12) # estimation basse
ts.plot(X,pred,pred_l,pred_u,xlab="t",ylab="Elec & Gas prod.", #sur graph complet
col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
ts.plot(window(X,start=c(2022,1)),pred,pred_l,pred_u,xlab="t", #sur dernière année
ylab="Elec & Gas prod.",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(1,3,2,2))
Nous affichons ici la série projetée en rouge et délimitons l’intervalle de confiance à 95% de cette série. L’aspect général de la série semble cohérent avec les données passées, mais il subsiste une incertitude quant à la précision de nos prévisions. De plus, il est notable que pour la période 2022-2023, la dernière observation de la série temporelle conserve une tendance similaire au cours des douze prochains instants.
Analyse a posteriori de la prévision
On tronque la série de l’année 2022, qu’on cherche ensuite à prévoir a partir de l’historique 1985-2021.
x_tronc=window(X,end=c(2021,12)) #tronquage de la série
y_tronc=log(x_tronc) # passage au log
x_a_prevoir=window(X,start=c(2022,1))On verifie que le modèle 5 sur la série tronquée est toujours validé.
model5tronc=Arima(y_tronc,order=c(1,1,1),list(order=c(1,1,2),period=12),include.mean=FALSE,method="CSS-ML") #vérification du modèle sur partie tronquée
summary(model5tronc)## Series: y_tronc
## ARIMA(1,1,1)(1,1,2)[12]
##
## Coefficients:
## ar1 ma1 sar1 sma1 sma2
## 0.4732 -0.9199 -0.4190 -0.2994 -0.4645
## s.e. 0.0528 0.0244 0.1983 0.1842 0.1483
##
## sigma^2 = 0.0006326: log likelihood = 970.75
## AIC=-1929.51 AICc=-1929.31 BIC=-1905.11
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.002087455 0.02463726 0.01938809 -0.04765542 0.4299135 0.6009786
## ACF1
## Training set -0.01456219On vérifie que les paramètres sont toujours significatifs.
test_model5tronc <- coeftest(model5) #test coeffs
test_model5tronc##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.465841 0.053352 8.7314 < 2.2e-16 ***
## ma1 -0.917634 0.025408 -36.1158 < 2.2e-16 ***
## sar1 -0.398903 0.189300 -2.1073 0.035095 *
## sma1 -0.326224 0.176746 -1.8457 0.064932 .
## sma2 -0.445371 0.143331 -3.1073 0.001888 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ils sont toujours aussi significatif que prédemment. On peut alors vérifier les résidus.
x5tronc<-model5tronc$residuals #test résidus
nlag=c(6,12,18,24,30,36)
res = c()
for (i in nlag){
a1=Box.test(x5tronc, lag=i, type = "Ljung-Box")
res=rbind(res,t(rbind(i,a1)[,3]))
}
colnames(res)=c("nlag","p-value x3")
res## nlag p-value x3
## [1,] 6 0.7401899
## [2,] 12 0.6534451
## [3,] 18 0.667341
## [4,] 24 0.4110009
## [5,] 30 0.385379
## [6,] 36 0.4745161shapiro.test(model5tronc$residuals)##
## Shapiro-Wilk normality test
##
## data: model5tronc$residuals
## W = 0.99502, p-value = 0.1641Les résidus peuveut donc être considérés comme un bruit blanc au niveau de test de 5% et ils vérifient aussi le test de normalité.
On peut maintenant superposer l’anée 2022 avec la prévision que nous avons obtenus en nous basant sur notre modèles des années 1985 à 2021.
pred_model5tronc=forecast(model5tronc,h=12,level=95) #prévision de l'année 2022
pred_tronc=exp(pred_model5tronc$mean)
pred_l_tronc=ts(exp(pred_model5tronc$lower),start=c(2022,1), #prévision basse
frequency=12)
pred_u_tronc=ts(exp(pred_model5tronc$upper),start=c(2022,1), # prévision haute
frequency=12)
ts.plot(x_a_prevoir,pred_tronc,pred_l_tronc,pred_u_tronc,xlab="t", ylab="Airpass",col=c(1,2,3,3),lty=c(1,1,2,2),lwd=c(3,3,2,2))
#supperposstion de la prévision et des données réels
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))
title(main = " Superposition de l’année 2022 avec les prévisions obtenus à partir de 1985 à 2021", line = 0.8,cex.main = 0.9)
On observe donc ici que la prévision en rouge est très proche de la réalisation en noir et que elle est bien contenu dans la l’intervalle de précision en vert.
Analyse des résultats
Notre étude initiale a révélé la présence de tendances et de saisonnalités dans la série temporelle, ce qui nous a incités à appliquer une différenciation d’ordre 2 pour rendre la série stationnaire.
Après avoir testé plusieurs modèles potentiels, nous avons finalement opté pour un modèle \(SARIMA(1,1,1)(1,1,1)_{12}\), qui a démontré des performances solides en tenant compte à la fois de la saisonnalité et de la tendance présentes dans les données.
L’estimation des paramètres du modèle choisi s’est avérée être un processus complexe, mais nous avons pu obtenir des valeurs cohérentes pour les coefficients AR, MA, SAR, et SMA, permettant ainsi de générer des prévisions significatives.
Nous avons validé la robustesse du modèle en vérifiant ses prévisions par rapport aux valeurs réelles tronquées de la série temporelle. Les critères d’erreur ont été utilisés pour évaluer la précision du modèle.
Conclusion
L’exercice s’est avéré être une expérience extrêmement instructive pour comprendre en profondeur la mise en œuvre des modèles de prévision ARMA, ARIMA, et SARIMA. La manipulation pratique des données réelles a permis de saisir les défis associés à la stationnarisation d’une série temporelle complexe et de sélectionner le modèle le mieux adapté à la série étudiée.
La prédiction précise des tendances saisonnières et des variations de la production électrique et de gaz aux États-Unis est d’une importance capitale pour la planification énergétique et économique. En se familiarisant avec les nuances et les difficultés réelles de l’application des modèles, nous avons acquis une compréhension plus profonde de leur pertinence et de leurs limites dans le contexte de la modélisation des séries temporelles complexes. Ainsi, pour interpreter nos résultats, nous pouvons voir qu’il y a une forte corrélation de production de gaz et electricité en hiver comme en été, nous pouvons surement lié celà entre autre à une forte consommation energétique pour le chauffage en période hivernale et à l’inverse, une baisse de cette consomation en été, 40% de l’énergie des foyers américains étant utilisée dans le chauffage d’après une étude de la US Energy Information Administration, disponible ici.
Cet exercice nous a également mis au défi de développer nos compétences en programmation en R et d’approfondir notre compréhension des concepts statistiques sous-jacents. En fin de compte, il a souligné l’importance de la rigueur analytique et de la validation adéquate des modèles pour assurer des prévisions fiables et précises.