Comment simuler à partir d'une copule gaussienne?

Supposons que j'ai deux distributions marginales univariées, disons $F$ et $G$ , à partir desquelles je peux simuler. Maintenant, construisons leur distribution conjointe en utilisant une copule gaussienne , notée $C(F,G;\Sigma)$ . Tous les paramètres sont connus.

Existe-t-il une méthode non-MCMC pour simuler à partir de cette copule?

Il existe une méthode très simple à simuler à partir de la copule gaussienne qui est basée sur les définitions de la distribution normale multivariée et de la copule de Gauss.

Je vais commencer par fournir la définition et les propriétés requises de la distribution normale multivariée, suivie de la copule gaussienne, puis je fournirai l'algorithme à simuler à partir de la copule de Gauss.

Distribution normale multivariée
Un vecteur aléatoire a une distribution normale multivariée si X d = μ + A Z , où Z est un vecteur k- dimensionnel de variables aléatoires normales normales indépendantes, μ est un d -vecteur dimensionnel de constantes, et A est une matrice d × k de constantes. La notation d$X = (X_1, \ldots, X_d)'$

$$X \stackrel{\mathrm{d}}{=} \mu + AZ,$$ $Z$$k$$\mu$$d$$A$$d\times k$$\stackrel{\mathrm{d}}{=}$dénote l'égalité dans la distribution. Ainsi, chaque composante de est essentiellement une somme pondérée de variables aléatoires normales standard indépendantes. A partir des propriétés des vecteurs moyens et des matrices de covariance, nous avons E ( X ) = μ et c o v ( X ) = Σ , avec Σ = A A ′ , conduisant à la notation naturelle X ∼ $X$
${\rm E}(X) = \mu$${\rm cov}(X) = \Sigma$$\Sigma = AA'$ .$X \sim {\rm N}_d(\mu, \Sigma)$

Copule de Gauss
La copule de Gauss est définie implicitement à partir de la distribution normale multivariée, c'est-à-dire que la copule de Gauss est la copule associée à une distribution normale multivariée. Plus précisément, d'après le théorème de Sklar, la copule de Gauss est où

$$C_P(u_1, \ldots, u_d) = \boldsymbol{\Phi}_P(\Phi^{-1}(u_1), \ldots, \Phi^{-1}(u_d)),$$ $\Phi$dénote la fonction de distribution normale standard, et désigne la fonction de distribution normale standard multivariée avec la matrice de corrélation P. Ainsi, la copule de Gauss est simplement une distribution normale multivariée standard où la transformation intégrale de probabilité est appliquée à chaque marge.$\boldsymbol{\Phi}_P$

Algorithme de simulation
Compte tenu de ce qui précède, une approche naturelle pour simuler à partir de la copule de Gauss consiste à simuler à partir de la distribution normale standard multivariée avec une matrice de corrélation appropriée , et à convertir chaque marge en utilisant la transformation intégrale de probabilité avec la fonction de distribution normale standard. Alors que la simulation à partir d'une distribution normale multivariée avec une matrice de covariance Σ revient essentiellement à faire une somme pondérée de variables aléatoires normales standard indépendantes, où la matrice de «poids» $P$$\Sigma$$A$ peut être obtenue par la décomposition de Cholesky de la matrice de covariance $\Sigma$ .

Par conséquent, un algorithme pour simuler échantillons de la copule de Gauss avec la matrice de corrélation P est:$n$$P$

1. Effectuer une décomposition Cholesky de et définir $P$$A$ comme matrice triangulaire inférieure résultante.
2. Répétez les étapes suivantes $n$ fois.
   1. Générer un vecteur $Z = (Z_1, \ldots, Z_d)'$ de variables normales standard indépendantes.
   2. Set $X = AZ$
   3. Retour .$U = (\Phi(X_1), \ldots, \Phi(X_d))'$

Le code suivant dans un exemple d'implémentation de cet algorithme utilisant R:

## Initialization and parameters 
set.seed(123)
P <- matrix(c(1, 0.1, 0.8,               # Correlation matrix
              0.1, 1, 0.4,
              0.8, 0.4, 1), nrow = 3)
d <- nrow(P)                             # Dimension
n <- 200                                 # Number of samples

## Simulation (non-vectorized version)
A <- t(chol(P))
U <- matrix(nrow = n, ncol = d)
for (i in 1:n){
    Z      <- rnorm(d)
    X      <- A%*%Z
    U[i, ] <- pnorm(X)
}

## Simulation (compact vectorized version) 
U <- pnorm(matrix(rnorm(n*d), ncol = d) %*% chol(P))

## Visualization
pairs(U, pch = 16,
      labels = sapply(1:d, function(i){as.expression(substitute(U[k], list(k = i)))}))

Le graphique suivant montre les données résultant du code R ci-dessus.