# clear environment variables
rm(list=ls())
# clear console
cat("\014")

##
library(readr)
library(zoo)
library(tseries)
library(forecast)
library(urca)
library(ggfortify)
library(caret)
library(stats)
library(lmtest)
library(car)
library(Metrics)
library(trend)
library(randtests)
library(uroot)

Auto Regressive Integrated Moving Average ARIMA

L’élévation du niveau de la mer est principalement causée par deux facteurs liés au réchauffement climatique : l’ajout d’eau provenant de la fonte des calottes glaciaires et des glaciers, et l’expansion de l’eau de mer à mesure qu’elle se réchauffe. Les données sur le niveau moyen global de la mer (GMSL) sont prises par les satellites de la NASA environ tous les 10 jours depuis 1993.

GMSL_variation_with_GIA : Variation du GMSL avec ajustement isostatique global (GIA)

L’ajustement isostatique global (GIA) ajuste les données pour refléter plus précisément les véritables changements du niveau de la mer, indépendamment des effets de ces mouvements terrestres.

Visualisation des données

data <- read_csv("data/GMSL_TPJAOS_5.1.csv")

## Rows: 1146 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (13): altimeter_type, merged_file_cycle, year_fraction, num_observations...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Générer une séquence de dates avec une fréquence de 10 jours, débutant en 1993
start_date <- as.Date("1993-01-01")
data$Date <- seq(start_date, by = "10 days", length.out = nrow(data))

# Trier les données par Date
data <- data[order(data$Date), ]

# Convertir en série temporelle (ts)
# Utiliser la colonne GMSL_variation_with_GIA
seaLevel <- ts(data$GMSL_variation_with_GIA, start=c(1993, 1), frequency=36.5)  # 36.5 périodes de 10 jours par an

# Afficher un aperçu de la série temporelle
plot(seaLevel)

## verifié si n'y a pas de valeurs manquantes
sum(is.na(seaLevel))

## [1] 0

Nous n’avons aucune valeurs manquantes dans la série

decomposeSeaLeavel <- decompose(seaLevel)
autoplot(decomposeSeaLeavel)

## Warning: Removed 70 rows containing missing values or values outside the scale range
## (`geom_line()`).

On soupçonne :

une tendance croissante :

Réalisation du test de Mann-Kindall pour vérification
```
mk.test(seaLevel)
```
```
## 
##  Mann-Kendall trend test
## 
## data:  seaLevel
## z = 45.269, n = 1146, p-value < 2.2e-16
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.857870e+05 1.674474e+08 8.929067e-01
```
Nous avons p-value < 5 % : Nous rejettons l’hypothèse nulle $H_0$ Cela indique une tendance significatif dans la série temporelles.

Effectuons un deuxième test afin d’être convaincu de la présence de tendance. Un test de Cox Stuart est effectué :
```
cox.stuart.test(seaLevel)
```
```
## 
##  Cox Stuart test
## 
## data:  seaLevel
## statistic = 573, n = 573, p-value < 2.2e-16
## alternative hypothesis: non randomness
```
p-value : < 5 % : Nous rejettons l’hypothèse nulle $H_0$ Cela indique une fois de plus une tendance monôtone significatif dans la série temporelles.

On conclut donc que la série comporte une tendance.

la présence de saisonnalité :

Nous effectuerons des tests pour vérification :

Effecutuons le test de Kruskall-Wallis.

kruskal.test(seaLevel~factor(cycle(seaLevel)))

## 
##  Kruskal-Wallis rank sum test
## 
## data:  seaLevel by factor(cycle(seaLevel))
## Kruskal-Wallis chi-squared = 20.047, df = 72, p-value = 1

p-value > 1 : Nous ne rejettons pas l’hypothèse null $H_0$. Nous avons une absence de saisonnalité dans notre série

oneway.test(seaLevel~factor(cycle(seaLevel)), var.equal = FALSE)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  seaLevel and factor(cycle(seaLevel))
## F = 0.22749, num df = 72.00, denom df = 367.26, p-value = 1

La p-value obtenue est extrêmement élevée (1), ce qui indique que nous ne rejetons pas l’hypothèse nulle. Cela signifie qu’il n’y a pas de différence significative entre les cycles, suggérant une absence de saisonnalité dans la série temporelle.

une certaine périodicité ou des motifs récurrents dans les résidus.

Nous ferons un test d’autocorralation des résidu de Box - Pierce
```
Box.test(decomposeSeaLeavel$random)
```
```
## 
##  Box-Pierce test
## 
## data:  decomposeSeaLeavel$random
## X-squared = 828.45, df = 1, p-value < 2.2e-16
```
Nous obtenons p-value < 5% : l’hypothèse nulle $H_0$ est ainsi rejeté. Nous avons donc une correlation des résidu.

Cela indique que bien que la tendance et la saisonnalité principales aient été retirées, il pourrait y avoir des composantes périodiques de plus haute fréquence qui n’ont pas été complètement capturées par le modèle de décomposition.

Nous pouvons donc conclure que la série comporte une tendance, une saison

Fonction d’autorrelation ACF

Acf(seaLevel)

Ici, on peut observer que la corrélation est très proche de 1 pour un décalage de 0 (corrélation parfaite avec elle-même). Ensuite, les valeurs d’autocorrélation décroissent lentement à mesure que le décalage augmente, mais restent significativement différentes de zéro, même pour des décalages importants.

Cela suggère une forte dépendance temporelle et une tendance potentielle dans les données de niveau de la mer. Une ACF décroissant lentement est souvent caractéristique des séries non stationnaires nécessitant un traitement supplémentaire (différenciation, etc.) avant une modélisation appropriée.

Fonction d’autorrelation partielle PACF

pacf(seaLevel)

On peut voir que la PACF est significative et proche de 1 pour un lag de 0, ce qui est normal puisqu’une série est parfaitement corrélée avec elle-même sans décalage.

Ensuite, la PACF diminue brusquement et semble devenir non significative (proche de zéro) pour les lags supérieurs à 1 ou 2. Cela suggère qu’après avoir pris en compte les dépendances aux lags proches, il n’y a plus de corrélation significative restante aux lags plus élevés.

Une PACF qui se coupe brusquement après quelques lags significatifs est souvent un indicateur d’un processus autorégressif d’ordre faible, comme un AR(1) ou AR(2), où les observations dépendent principalement des quelques observations précédentes.

Verification de la non stationnarité de la série

adf.test(seaLevel, k=36.5)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  seaLevel
## Dickey-Fuller = -1.5779, Lag order = 36.5, p-value = 0.757
## alternative hypothesis: stationary

p-value > 0.05 : On ne peut rejetter ainsi l’hypothèse nulle. Donc la série est belle et bien non stationnaire. On doit donc différencier la série pour la rendre stationnaire.

Identification des paramètres du modèle

Ordre de differentiation (d)

d = ndiffs(seaLevel)
d

## [1] 1

Il faut donc différencier la série une seule fois pour la rendre stationnaire. Donc notre nouvel dataset est :

seaLevelDiff = diff(seaLevel)

Ordre (p)

pacf(seaLevelDiff)

D’après le graphe de la fonction d’autocorrélation partielle (PACF) pour la série “seaLevelDiff”, on peut suggérer un ordre autorégressif p=1 ou p=2 ou p=3 pour un éventuel modèle ARIMA.

Ordre (q)

Acf(seaLevelDiff)

Les 2 premiers lags semblent avoir des valeurs sigificatifs (sortant des bandes de confiance). On peut prendre q=1 ou q=2.

Préparation des Données

n = length(seaLevel)
train_data = seaLevel[1:round(0.8*n)]
test_data = seaLevel[(round(0.8*n) + 1):1126]

On divise la serie seaLevel en ensembles d’entraînement et de test selon un ratio de 80/20, où 80% des données sont utilisées pour l’entraînement et 20% pour les tests.

Modèle

best_model = NULL
param_p = 0
param_q = 0
aic_best_criterion_info = 1000000
bic_best_criterion_info = 1000000

for (p in 1:3) {
  for (q in 1:3) {
    # Créer le modèle ARIMA
    
    model <- Arima(train_data, order = c(p, 1, q), include.drift=TRUE)
    aic_criterion_info = AIC(model)
    bic_criterion_info = BIC(model)
    cat("ARIMA (", p, ', 1 ,',  q, ') : AIC : ', aic_criterion_info, "BIC : ", bic_criterion_info ,  "\n")
    if (aic_criterion_info < aic_best_criterion_info && bic_criterion_info < bic_best_criterion_info){
     best_model = model
     aic_best_criterion_info = aic_criterion_info
     bic_best_criterion_info = bic_criterion_info
     param_q = q
     param_p = p
    }
  }
}

## ARIMA ( 1 , 1 , 1 ) : AIC :  4257.12 BIC :  4276.4 
## ARIMA ( 1 , 1 , 2 ) : AIC :  4252.129 BIC :  4276.229 
## ARIMA ( 1 , 1 , 3 ) : AIC :  4252.231 BIC :  4281.151 
## ARIMA ( 2 , 1 , 1 ) : AIC :  4249.301 BIC :  4273.401 
## ARIMA ( 2 , 1 , 2 ) : AIC :  4224.236 BIC :  4253.156 
## ARIMA ( 2 , 1 , 3 ) : AIC :  4213.754 BIC :  4247.494 
## ARIMA ( 3 , 1 , 1 ) : AIC :  4248.859 BIC :  4277.779 
## ARIMA ( 3 , 1 , 2 ) : AIC :  4217.232 BIC :  4250.972 
## ARIMA ( 3 , 1 , 3 ) : AIC :  4226.782 BIC :  4265.342

cat("Le meilleur modèle est obtenu est un ARIMA (", param_p, ', 1 ,', param_q, ')')

## Le meilleur modèle est obtenu est un ARIMA ( 2 , 1 , 3 )

Le meilleur modèle minimissant les criteres AIC et BIC est celui avec p=2 et q=3 ainsi on a un ARIMA(2, 1, 3).

Determination des coefficients du modèle

summary(best_model)

## Series: train_data 
## ARIMA(2,1,3) with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2      ma3   drift
##       0.9718  -0.9912  -1.1534  1.1345  -0.1556  0.0978
## s.e.  0.0081   0.0077   0.0376  0.0443   0.0365  0.0640
## 
## sigma^2 = 5.763:  log likelihood = -2099.88
## AIC=4213.75   AICc=4213.88   BIC=4247.49
## 
## Training set error measures:
##                        ME    RMSE      MAE  MPE MAPE     MASE       ACF1
## Training set -0.001897246 2.39137 1.841836 -Inf  Inf 0.958979 0.01375843

Significativité des coefficients :

Pour déterminer la significativité des coefficients, nous calculons les ratios (coefficient / erreur standard) et comparons ces valeurs aux seuils critiques de la distribution normale.

Calcul des ratios :

ar1 : 0.9706 / 0.0083 $\approx$ 119,97 > 1.96
ar2 : -0.9912 / 0.0077 $\approx$ 128,72 > 1.96
ma1 : -1.1534 / 0.0443 $\approx$ -30,675 < -1.96
ma2 : 1.1345 / 0.0443 $\approx$ 6.17 > 25,609
ma3 : -0.1556 / 0.0365 $\approx$ -4,263 < -1.96

Tous les coefficients (ar1, ar2, ma1, ma2 et ma3) sont significatifs car leurs ratios (valeurs t) sont tous supérieurs en valeur absolue à 1.96. Cela signifie que ces coefficients sont significativement différents de zéro à un niveau de confiance de 95%.

Validation du Modèle

Affichage des résidus :

residuals = best_model$residuals
autoplot(residuals)

Stationnarité des résidus

adf.test(residuals, k=36.5)

## Warning in adf.test(residuals, k = 36.5): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  residuals
## Dickey-Fuller = -7.1407, Lag order = 36.5, p-value = 0.01
## alternative hypothesis: stationary

On a p-value < 0.05 alors on rejette l’hypothèse nulle. Les résidus sont stationnaires.

Test d’autocorrélation

plot(resid(best_model), type = "l")

Test de Box et Pierce

Box.test(residuals)

## 
##  Box-Pierce test
## 
## data:  residuals
## X-squared = 0.17358, df = 1, p-value = 0.6769

Le résultat du test de Box-Pierce sur les résidus montre une valeur de la statistique du chi-deux (X-squared) de 0.17358 avec 1 degré de liberté (df). La p-value associée est de 0,6769. Avec une p-value bien supérieure au seuil habituel de significativité de 0.05, on ne peut pas rejeter l’hypothèse nulle d’absence d’autocorrélation dans les résidus. En d’autres termes, ce test ne détecte pas de dépendance significative dans les résidus du modèle, ce qui est généralement un bon signe. Cela suggère que le modèle a bien capturé les dépendances temporelles dans les données et qu’il ne reste pas d’autocorrélation exploitable dans les résidus.

Test de normalité des résidus

hist(residuals)

On soupçonne la normalité des résidus. On effectue un test de normalité pour vérifier.

shapiro.test(residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals
## W = 0.98911, p-value = 2.551e-06

On a p-value < 0.05. On rejette alors l’hypothese nulle. On conclut alors que les résidus ne suivent pas une distribution normale.

Prédiction

Après avoir ajusté un modèle ARIMA à nos données et évalué la significativité de ses coefficients, nous sommes maintenant prêts à utiliser ce modèle pour faire des prédictions.

forecastSeaLevel = forecast(best_model, level=c(95), h=72)
autoplot(forecastSeaLevel)

Le modèle capture bien la tendance à la hausse des données historiques.
L’incertitude augmente au fur et à mesure que l’on s’éloigne dans le temps, ce qui est typique des prévisions de séries temporelles. Plus l’horizon de prévision est lointain, plus l’intervalle de confiance est large, indiquant une plus grande incertitude.

Nous utiliserons, la fonction auto.arima de la bibliothèque forecast en R pour identifier les meilleurs paramètres du modèle ARIMA pour la série temporelle .

autoArimaModel = auto.arima(train_data)
summary(autoArimaModel)

## Series: train_data 
## ARIMA(2,1,2) with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2   drift
##       0.9408  -0.9069  -1.0307  0.8884  0.0983
## s.e.  0.0267   0.0403   0.0256  0.0400  0.0707
## 
## sigma^2 = 5.844:  log likelihood = -2106.12
## AIC=4224.24   AICc=4224.33   BIC=4253.16
## 
## Training set error measures:
##                        ME     RMSE     MAE  MPE MAPE      MASE        ACF1
## Training set -0.002036341 2.409576 1.85285 -Inf  Inf 0.9647136 -0.06571887

Les paramètres $p$ , $d$ et $q$ obtenus par la fonction auto.arima semblent presque identique aux paramètres obtenu à partir de l’ACF et du PACF.

Vérification de l’autocorrelation des erreurs

plot(resid(autoArimaModel), type = "l")

Nous utiliserons un test de Box-Pierce pour cette vérification

Box.test(resid(autoArimaModel))

## 
##  Box-Pierce test
## 
## data:  resid(autoArimaModel)
## X-squared = 3.9605, df = 1, p-value = 0.04658

Nous avons une p-value sensiblement égale à 5%. Nous effectuerons un autre test (test de Durbin Wason) pour y voir claire

dwtest(resid(autoArimaModel) ~ 1)

## 
##  Durbin-Watson test
## 
## data:  resid(autoArimaModel) ~ 1
## DW = 2.1309, p-value = 0.9764
## alternative hypothesis: true autocorrelation is greater than 0

La p-value associée au test est très élevée (0.9764), bien au-dessus du seuil de 0.05 généralement utilisé pour la signification statistique. Une p-value élevée indique que nous ne pouvons pas rejeter l’hypothèse nulle , qui postule l’absence d’autocorrélation des résidus.

Par conséquent, on peut conclure que les résidus de votre modèle ARIMA se comportent de manière indépendante, ce qui est un bon signe de la qualité du modèle ajusté.

Normalité des résidus

shapiro.test(autoArimaModel$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  autoArimaModel$residuals
## W = 0.98965, p-value = 4.631e-06