par=(mfrow=c(1,3))

plot(df_var$Rainfall_mm, df_var$SoilMoisture_pct,
     xlab = "Précipitations (mm)",
     ylab = "Humidité du sol (%)",
     main = "Nuage de points pluie – humidité du sol")

plot(df_var$Rainfall_mm, df_var$WaterLevel_m,
     xlab = "Précipitations (mm)",
     ylab = "Niveau de l'eau (m)",
     main = "Nuage de points pluie – Niveau de l'eau")

plot(df_var$WaterLevel_m, df_var$SoilMoisture_pct,
     xlab = "Niveau de l'eau (m)",
     ylab = "Humidité du sol (%)",
     main = "Niveau de l'eau – humidité du sol")

• Sur le premier graphique, l’on observe une dépendance monotone entre les précipitations et le taux d’humidité du sol avec une forte dépendance dans la queue de distribution en haut à droite

• Sur le second graphique, l’on constate que les faibles niveaux de pluie (-20mm) n’ont pas d’effet systématique sur le niveau de cours d’eau. Mais au delà de 20 mm de pluie, l’on constate une dépendance monotone entre les précipitations et le niveau du cours d’eau avec une forte dépendance dans la queue en haut à droite

• Sur le troisième graphique, l’on ne constate ne soupçonne pas de lien apparent dépendance entre les deux variables. Nous sommes proche de la quasi-indépendance

  Dans la suite du travail, nous prenons la peine de vérifier ces observations avec des approches beaucoup plus rigoureuses vues en cours

cor.test(x=U_mat1[,1], y=U_mat1[,2], method="kendall")

## 
##  Kendall's rank correlation tau
## 
## data:  U_mat1[, 1] and U_mat1[, 2]
## z = 21.915, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.3440394

cor.test(x=U_mat2[,1], y=U_mat2[,2], method="kendall")

## 
##  Kendall's rank correlation tau
## 
## data:  U_mat2[, 1] and U_mat2[, 2]
## z = 21.702, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.3585618

cor.test(x=U_mat3[,1], y=U_mat3[,2], method="kendall")

## 
##  Kendall's rank correlation tau
## 
## data:  U_mat3[, 1] and U_mat3[, 2]
## z = 12.112, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.2002208

D’après les tests, il existe une relation de dépendance entre les deux variables de chaque paire. L’on remarque que ce sont des dépendances relativement faibles au regards des taux de Kendall respectifs estimés.

En particulier, pour le troisième couple, contrairement à notre intuition, l’on a bel et bien une relation de dépendance faible mais significative ! Mais le problème est que la forme du nuage de points ne nous permet pas d’identifier aisément le type de copule à laquelle l’on a affaire! On verra ça.

Le test de Kendall nous permet d’apprècier l’existence d’une relation de dépendance générale entre les variables. Mais nous avons aussi besoin de confirmer nos observations au niveau des queues. Pour cela, nous allons réaliser un x-plot ( lire chi-plot ) à partir des pseudo-observations. En plus, le chi-plot peut permettre aussi de trancher au niveau de la dépendance globale.

#install.packages("VineCopula")
library(VineCopula)

BiCopChiPlot(u1 = U_mat1[,1], u2 = U_mat1[,2], PLOT = T)

• Les valeurs de chi sont toutes positives. Cela confirme une dépendance globale positive.

• Au centre ( autour de lamda=0), l’on a des valeurs de chi élevées. l’on a donc une dépendance forte au centre. C’est-à-dire les valeurs typiques tendent à évoluer de manière concordantes

• L’on remarque également une dépendance forte au niveau supérieur droit car pour lambda proche de 1, l’on a une forte concentration et l’on a les plus hautes valeurs de chi

Toutes ces conditions nous font penser à une copule de Gumbel ou une copule de survie de Clayton

BiCopChiPlot(u1 = U_mat2[,1], u2 = U_mat2[,2], PLOT = T)

• De même l’on a une dépendance positive mais faible en générale car les valeurs de chi bien qu’étant positive, sont comprises, pour la plupart entre 0 et 0.5 avec une dépendance relativement forte au centre.

• Elle semble laisser paraître une forte dépendance de queue à droite

  Ce graphe a une structure assez capricieuse. Nous sommes réticents mais au regard des pseudo observations, nous allons miser sur une copule de Gumbel ou une copule de survie de Clayton

BiCopChiPlot(u1 = U_mat3[,1], u2 = U_mat3[,2], PLOT = T)

De façon générale nous avons des valeurs de chi positive et alignées en bande le long de l’axe chi=0. Alors la dépendance est positive mais elle est trop faible. De plus, en raison de la structure horizontale le long de l’axe des abscisses, l’on se rapproche ici d’une situation d’indépendance . Donc nous n’allons pas estimer de copules pour ce couple de variable dans la suite de notre travail.

# Coefficient de queue de droite
lambda_U_hat <- function(U, q ) {
  u <- U[,1]
  v <- U[,2]
  num <- sum(u > 1-q & v > 1-q)
  den <- sum(u > 1-q)
  if (den == 0) return(NA_real_)
  num / den
}

# coefficient de queue de gauche
lambda_L_hat <- function(U, q ) {
  u <- U[,1]
  v <- U[,2]
  num <- sum(u <= q & v <= q)
  den <- sum(u <= q)
  if (den == 0) return(NA_real_)
  num / den
}

q_grid <- seq(0.01, 0.99, by = 0.01)

par=(mfrow=c(1,2))
# estimation en queue supérieure

lambda_val1 <- sapply(q_grid, function(q) lambda_U_hat(U_mat1, q))
plot(q_grid, lambda_val1,
     type = "l",
     lwd = 2,
     xlab = "q",
     ylab = expression(hat(lambda)[U](q)),
     main = "Estimateur empirique de la dépendance en queue supérieure")

abline(h = 0, lty = 2, col = "grey")

# Estimateur en queue inférieure

lambda_val2 <- sapply(q_grid, function(q) lambda_L_hat(U_mat1, q))
plot(q_grid, lambda_val2,
     type = "l",
     lwd = 2,
     xlab = "q",
     ylab = expression(hat(lambda)[L](q)),
     main = "Estimateur empirique de la dépendance en queue inférieure")

abline(h = 0, lty = 2, col = "grey")

• à gauche λL = 0

• à droite λU = 1

par=(mfrow=c(1,2))
# estimation en queue supérieure

lambda_val3 <- sapply(q_grid, function(q) lambda_U_hat(U_mat2, q))
plot(q_grid, lambda_val3,
     type = "l",
     lwd = 2,
     xlab = "q",
     ylab = expression(hat(lambda)[U](q)),
     main = "Estimateur empirique de la dépendance en queue supérieure")

abline(h = 0, lty = 2, col = "grey")

# Estimateur en queue inférieure

lambda_val4 <- sapply(q_grid, function(q) lambda_L_hat(U_mat2, q))
plot(q_grid, lambda_val4,
     type = "l",
     lwd = 2,
     xlab = "q",
     ylab = expression(hat(lambda)[L](q)),
     main = "Estimateur empirique de la dépendance en queue inférieure")

abline(h = 0, lty = 2, col = "grey")

• à gauche λL= 0

• à droite λU = 1

Nous allons à présent estimer le paramètre de dépendance θ pour chaque couple. En fait, le paramètre de dépendance contrôle la force et la structure de dépendance entre les variables. Pour faire simple, il détermine ” à quel point les variables sont liées”

#estimation de tau pour Gumbel
tau_emp1 <- cor(U_mat1[,1], U_mat1[,2], method = "kendall") 

#estimation de tau pour la copule de survie de Clayton
U_surv1 <- 1 - U_mat1
tau_surv1 <- cor(U_surv1[,1], U_surv1[,2], method = "kendall") 

theta_gumb_mm1 <- iTau(gumbelCopula(), tau_emp1)
theta_sclay_mm1 <- iTau(claytonCopula(), tau_surv1)  

theta_gumb_mm1; theta_sclay_mm1

## [1] 1.524482

## [1] 1.048963

tau_emp2 <- cor(U_mat2[,1], U_mat2[,2], method = "kendall") #estimation de tau pour Gumbel

U_surv2 <- 1 - U_mat2
tau_surv2 <- cor(U_surv2[,1], U_surv2[,2], method = "kendall") #estimation de tau pour la copule de survie de Clayton

theta_gumb_mm2 <- iTau(gumbelCopula(), tau_emp2)
theta_sclay_mm2 <- iTau(claytonCopula(), tau_surv2)  

theta_gumb_mm2; theta_sclay_mm2

## [1] 1.558997

## [1] 1.117993

# Gumbel
fitGumb1 <- fitCopula(
  gumbelCopula(dim = 2),
  U_mat1,
  method = "ml"
)

# Clayton de survie
U_surv1 <- 1 - U_mat1
fitSClay1 <- fitCopula(
  claytonCopula(dim = 2),
  U_surv1,
  method = "ml"
)

summary(fitGumb1)

## Call: fitCopula(gumbelCopula(dim = 2), U_mat1, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Gumbel copula, dim. d = 2 
##       Estimate Std. Error
## alpha    1.486      0.027
## The maximized loglikelihood is 284.1 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        6        6

summary(fitSClay1)

## Call: fitCopula(claytonCopula(dim = 2), U_surv1, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Clayton copula, dim. d = 2 
##       Estimate Std. Error
## alpha    1.049      0.046
## The maximized loglikelihood is 266.3 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        3        3

# Gumbel
fitGumb2 <- fitCopula(
  gumbelCopula(dim = 2),
  U_mat2,
  method = "ml"
)

# Clayton de survie
U_surv2 <- 1 - U_mat2
fitSClay2 <- fitCopula(
  claytonCopula(dim = 2),
  U_surv2,
  method = "ml"
)

summary(fitGumb2)

## Call: fitCopula(gumbelCopula(dim = 2), U_mat2, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Gumbel copula, dim. d = 2 
##       Estimate Std. Error
## alpha    1.494      0.028
## The maximized loglikelihood is 277.1 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        7        7

summary(fitSClay2)

## Call: fitCopula(claytonCopula(dim = 2), U_surv2, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Clayton copula, dim. d = 2 
##       Estimate Std. Error
## alpha    1.118      0.047
## The maximized loglikelihood is 251.1 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        3        3

Pour chacun des couples, les deux méthodes donnent des estimations proches pour les deux copules. On voit que pour les deux copules la valeur de l’estimée de θ montre bien qu’il y a une dépendence.

Nous allons maintenant comparer les deux modèles à l’aide de l’AIC/BIC (possible pour l’estimation par maximum de vraisemblance)

logLik(fitGumb1); logLik(fitSClay1)

## 'log Lik.' 284.0747 (df=1)

## 'log Lik.' 266.3464 (df=1)

logLik(fitGumb2); logLik(fitSClay2)

## 'log Lik.' 277.0811 (df=1)

## 'log Lik.' 251.132 (df=1)

Il ressort que la copule de Gumbel est le meilleur modèle dans les deux cas

Test de Cramér-Von Mises

Comparaison des lignes de niveau:

Copule de Gumbel

library(lattice)
library(latticeExtra)


# Calculs
fitGumb1 <- fitCopula(gumbelCopula(), as.matrix(U_mat1), method = "mpl")
emp_cop1 <- empCopula(as.matrix(U_mat1))

# Graphique
contour(fitGumb1@copula, pCopula, 
        main = "Gumbel : Théorique (noir) vs Empirique (rouge)",
        col = "black")

# Superposition
contour(emp_cop1, pCopula, 
        col = "red", 
        lty = 2, 
        add = TRUE)

# Légende (Correction du symbole inattendu ici)
legend("bottomright", legend = c("Théorique", "Empirique"), 
       col = c("black", "red"), lty = c(1, 2), lwd = 1.5)

Copule de Clayton

# 1. Préparation de la copule théorique (Clayton)
# On récupère l'objet copule ajusté
clay_cop <- fitSClay1@copula

# 2. Création du graphique de base avec la copule théorique
# On utilise contour() qui est la méthode native et stable
contour(clay_cop, pCopula, 
        main = "Clayton : Théorique (noir) vs Empirique (rouge)",
        col = "black", 
        nlevels = 20)

# 3. Calcul de la copule empirique sur une grille
u <- seq(0, 1, length.out = 30) # Plus de points pour un contour plus lisse
grid <- as.matrix(expand.grid(u1 = u, u2 = u))

# Calcul des valeurs C.n (Copule empirique)
# C.n est la fonction de répartition empirique du package copula
val_emp <- C.n(grid, X = as.matrix(U_mat1))

# Transformation en matrice pour la fonction contour() de base
z_matrix <- matrix(val_emp, nrow = length(u))

# 4. Superposition de la copule empirique
# add = TRUE permet de dessiner par-dessus le premier graphique
contour(u, u, z_matrix, 
        col = "red", 
        lty = 2, 
        add = TRUE, 
        nlevels = 20, 
        drawlabels = FALSE)

# 5. Légende
legend("bottomright", legend = c("Clayton Théorique", "Empirique"),
       col = c("black", "red"), lty = c(1, 2), bty = "n")

Test de Cramér-Von Mises

Comparaison des lignes de niveau:

Copule de Gumbel

emp_cop2 <- empCopula(as.matrix(U_mat2))

# 1. Copule théorique
contour(fitGumb2@copula, pCopula, 
        main = "Gumbel : Théorique (noir) vs Empirique (rouge)",
        col = "black", 
        lwd = 1.5)

# 2. Copule empirique (superposition)
contour(emp_cop2, pCopula, 
        col = "red", 
        lty = 2, 
        add = TRUE)

# 3. Légende
legend("bottomright",
       legend = c("Théorique", "Empirique"),
       col = c("black", "red"),
       lty = c(1, 2),
       lwd = 1.5)

Copule de Clayton

# 1. Préparation de la copule théorique (Clayton)
# On récupère l'objet copule ajusté
clay_cop <- fitSClay2@copula

# 2. Création du graphique de base avec la copule théorique
# On utilise contour() qui est la méthode native et stable
contour(clay_cop, pCopula, 
        main = "Clayton : Théorique (noir) vs Empirique (rouge)",
        col = "black", 
        nlevels = 20)

# 3. Calcul de la copule empirique sur une grille
u <- seq(0, 1, length.out = 30) # Plus de points pour un contour plus lisse
grid <- as.matrix(expand.grid(u1 = u, u2 = u))

# Calcul des valeurs C.n (Copule empirique)
# C.n est la fonction de répartition empirique du package copula
val_emp <- C.n(grid, X = as.matrix(U_mat2))

# Transformation en matrice pour la fonction contour() de base
z_matrix <- matrix(val_emp, nrow = length(u))

# 4. Superposition de la copule empirique
# add = TRUE permet de dessiner par-dessus le premier graphique
contour(u, u, z_matrix, 
        col = "red", 
        lty = 2, 
        add = TRUE, 
        nlevels = 20, 
        drawlabels = FALSE)

# 5. Légende
legend("bottomright", legend = c("Clayton Théorique", "Empirique"),
       col = c("black", "red"), lty = c(1, 2), bty = "n")