¿Cómo simular desde una cópula gaussiana?

Supongamos que tengo dos distribuciones marginales univariadas, digamos $F$ y $G$ , que puedo simular. Ahora, construya su distribución conjunta usando una cópula gaussiana , denotada $C(F,G;\Sigma)$ . Todos los parámetros son conocidos.

¿Existe algún método que no sea MCMC para simular a partir de esta cópula?

Existe un método muy simple para simular a partir de la cópula gaussiana que se basa en las definiciones de la distribución normal multivariada y la cópula Gauss.

Comenzaré proporcionando la definición y las propiedades requeridas de la distribución normal multivariada, seguida de la cópula gaussiana, y luego proporcionaré el algoritmo para simular a partir de la cópula Gauss.

Distribución normal multivariante
Un vector aleatorio tiene una distribución normal multivariante si X d = μ + A Z , donde Z es un vector k- dimensional de variables aleatorias normales estándar independientes, μ es un vector d- dimensional de constantes, y A es una matriz de constantes d × k . La notación d$X = (X_1, \ldots, X_d)'$

$$X \stackrel{\mathrm{d}}{=} \mu + AZ,$$ $Z$$k$$\mu$$d$$A$$d\times k$$\stackrel{\mathrm{d}}{=}$denota igualdad en la distribución. Por lo tanto, cada componente de es esencialmente una suma ponderada de variables aleatorias normales estándar independientes. A partir de las propiedades de los vectores medios y las matrices de covarianza, tenemos E ( X ) = μ y c o v ( X ) = Σ , con Σ = A A ' , lo que lleva a la notación natural X ∼ N d ( μ , Σ ) .$X$
${\rm E}(X) = \mu$${\rm cov}(X) = \Sigma$$\Sigma = AA'$$X \sim {\rm N}_d(\mu, \Sigma)$

Cópula de Gauss
La cópula de Gauss se define implícitamente a partir de la distribución normal multivariada, es decir, la cópula de Gauss es la cópula asociada con una distribución normal multivariada. Específicamente, a partir del teorema de Sklar, la cópula de Gauss es donde

$$C_P(u_1, \ldots, u_d) = \boldsymbol{\Phi}_P(\Phi^{-1}(u_1), \ldots, \Phi^{-1}(u_d)),$$ $\Phi$denota la función de distribución normal estándar, y denota la función de distribución normal estándar multivariada con la matriz de correlación P. Entonces, la cópula de Gauss es simplemente una distribución normal multivariada estándar donde la transformación integral de probabilidad se aplica a cada margen.$\boldsymbol{\Phi}_P$

Algoritmo de simulación
En vista de lo anterior, un enfoque natural para simular desde la cópula de Gauss es simular desde la distribución normal estándar multivariada con una matriz de correlación apropiada , y convertir cada margen usando la transformación integral de probabilidad con la función de distribución normal estándar. Mientras que la simulación de una distribución normal multivariada con matriz de covarianza Σ esencialmente se reduce a hacer una suma ponderada de variables aleatorias normales estándar independientes, donde la matriz de "peso" A se puede obtener mediante la descomposición de Cholesky de la matriz de covarianza $P$$\Sigma$$A$$\Sigma$ .

Por lo tanto, un algoritmo para simular muestras de la cópula de Gauss con la matriz de correlación P es:$n$$P$

1. Realice una descomposición Cholesky de y establezca $P$$A$ como la matriz triangular inferior resultante.
2. Repita los siguientes pasos $n$ veces.
   1. Genere un vector $Z = (Z_1, \ldots, Z_d)'$ de variables normales estándar independientes.
   2. Conjunto $X = AZ$
   3. Devuelve .$U = (\Phi(X_1), \ldots, \Phi(X_d))'$

El siguiente código en una implementación de ejemplo de este algoritmo usando 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)))}))

El siguiente cuadro muestra los datos resultantes del código R anterior.