Como simular a partir de uma cópula gaussiana?

Suponha que eu tenha duas distribuições marginais univariadas, digamos e , das quais eu possa simular. Agora, construa sua distribuição conjunta usando uma cópula gaussiana , denominada . Todos os parâmetros são conhecidos.$F$$G$$C(F,G;\Sigma)$

Existe um método não MCMC para simular a partir desta cópula?

Existe um método muito simples de simular a partir da cópula gaussiana, que se baseia nas definições da distribuição normal multivariada e na cópula de Gauss.

Começarei fornecendo a definição e as propriedades necessárias da distribuição normal multivariada, seguida pela cópula gaussiana, e depois fornecerei o algoritmo para simular a partir da cópula de Gauss.

Distribuição normal multivariada
Um vetor aleatório tem uma distribuição normal multivariada se X d = μ + A Z , onde Z é um vetor tridimensional k de variáveis ​​aleatórias normais padrão independentes, μ é um vetor d- dimensional de constantes e A é uma matriz d × k de constantes. A notação 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 igualdade na distribuição. Portanto, cada componente de é essencialmente uma soma ponderada de variáveis ​​aleatórias normais padrão independentes. A partir das propriedades dos vetores médios e matrizes de covariância, temos E ( X ) = μ e c o v ( X ) = Σ , com Σ = A A ′ , levando à notação 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
A cópula de Gauss é definida implicitamente a partir da distribuição normal multivariada, ou seja, a cópula de Gauss é a cópula associada a uma distribuição normal multivariada. Especificamente, a partir do teorema de Sklar, a cópula de Gauss é onde

$$C_P(u_1, \ldots, u_d) = \boldsymbol{\Phi}_P(\Phi^{-1}(u_1), \ldots, \Phi^{-1}(u_d)),$$ $\Phi$denota a função de distribuição normal padrão e denota a função de distribuição normal padrão multivariada com matriz de correlação P. Portanto, a cópula de Gauss é simplesmente uma distribuição normal multivariada padrão em que a transformação integral de probabilidade é aplicada a cada margem.$\boldsymbol{\Phi}_P$

Algoritmo de simulação
Em vista do exposto, uma abordagem natural para simular a partir da cópula de Gauss é simular a partir da distribuição normal padrão multivariada com uma matriz de correlação apropriada e converter cada margem usando a transformação integral de probabilidade com a função de distribuição normal padrão. Embora simulando a partir de uma distribuição normal multivariada com matriz covariância Σ essencialmente vem para baixo para fazer uma soma ponderada de variáveis independentes normais padrão aleatório, em que o "peso" matriz A podem ser obtidos por decomposição de Cholesky da matriz covariância Σ .$P$$\Sigma$$A$$\Sigma$

Portanto, um algoritmo para simular amostras da cópula de Gauss com matriz de correlação P é:$n$$P$

1. Realize uma decomposição de de Cholesky e defina A como a matriz triangular inferior resultante.$P$$A$
2. Repita as etapas a seguir vezes. $n$
   1. Gere um vetor de variáveis ​​normais padrão independentes.$Z = (Z_1, \ldots, Z_d)'$
   2. Defina $X = AZ$
   3. Retorno .$U = (\Phi(X_1), \ldots, \Phi(X_d))'$

O código a seguir em um exemplo de implementação desse 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)))}))

A tabela a seguir mostra os dados resultantes do código R acima.