Wie simuliert man eine Gaußsche Kopula?

Angenommen, ich habe zwei univariate Randverteilungen, sagen wir und , die ich simulieren kann. Konstruieren Sie nun ihre gemeinsame Verteilung mit einer Gaußschen Kopula , die mit . Alle Parameter sind bekannt.G C ( F , G ; Σ $F$$G$$C(F,G;\Sigma)$

Gibt es eine Nicht-MCMC-Methode zum Simulieren von dieser Kopula?

Es gibt eine sehr einfache Methode, um aus der Gaußschen Copula zu simulieren, die auf den Definitionen der multivariaten Normalverteilung und der Gaußschen Copula basiert.

Zunächst werde ich die erforderliche Definition und die Eigenschaften der multivariaten Normalverteilung, gefolgt von der Gaußschen Copula, bereitstellen und dann den Algorithmus zum Simulieren aus der Gaußschen Copula bereitstellen.

Multivariate Normalverteilung
Ein Zufallsvektor hat eine multivariate Normalverteilung, wenn X d = μ + A Z ist , wobei Z ein k- dimensionaler Vektor unabhängiger normaler Zufallsvariablen ist, μ a ist d- dimensionaler Konstantenvektor, und A ist eine d × k- Konstantenmatrix. Die 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}}{=}$bezeichnet Gleichheit in der Verteilung. Jede Komponente von ist also im Wesentlichen eine gewichtete Summe unabhängiger normaler Standard-Zufallsvariablen. Aus den Eigenschaften der mittleren Vektoren und Kovarianzmatrizen ergibt sich E ( X ) = μ und c o v ( X ) = Σ , wobei Σ = A A ′ , was zur natürlichen Notation X ∼ N d ( μ , Σ ) führt .$X$
${\rm E}(X) = \mu$${\rm cov}(X) = \Sigma$$\Sigma = AA'$$X \sim {\rm N}_d(\mu, \Sigma)$

Gauß-Copula
Die Gauß-Copula wird implizit aus der multivariaten Normalverteilung definiert, dh die Gauß-Copula ist die Copula, die einer multivariaten Normalverteilung zugeordnet ist. Insbesondere aus Theorem Sklar ist die Gauss - Copula ist wobei die Standardnormalverteilungsfunktion bezeichnet und die multivariate Standardnormalverteilungsfunktion mit der Korrelationsmatrix P bezeichnet. Die Gauß-Kopula ist also einfach eine multivariate Standardnormalverteilung mit der Wahrscheinlichkeitsintegraltransformation wird auf jeden Rand angewendet.Φ Φ

$$C_P(u_1, \ldots, u_d) = \boldsymbol{\Phi}_P(\Phi^{-1}(u_1), \ldots, \Phi^{-1}(u_d)),$$ $\Phi$$\boldsymbol{\Phi}_P$

Simulationsalgorithmus
Angesichts des Vorstehenden besteht ein natürlicher Ansatz zum Simulieren aus der Gauß-Kopula darin, aus der multivariaten Standardnormalverteilung mit einer geeigneten Korrelationsmatrix zu simulieren und jeden Rand unter Verwendung der Wahrscheinlichkeitsintegraltransformation mit der Standardnormalverteilungsfunktion umzuwandeln. Bei der Simulation einer multivariaten Normalverteilung mit der Kovarianzmatrix Wesentlichen eine gewichtete Summe unabhängiger normaler Zufallsvariablen berechnet, wobei die "Gewichts" -Matrix durch die Cholesky-Zerlegung der Kovarianzmatrix .Σ A $P$$\Sigma$$A$$\Sigma$

Daher ist ein Algorithmus zum Simulieren von Abtastwerten aus der Gauß-Kopula mit der Korrelationsmatrix :$n$$P$

1. Führen Sie eine Cholesky-Zerlegung von durch und setzen Sie als resultierende untere Dreiecksmatrix.$P$$A$
2. Wiederholen Sie die folgenden Schritte Mal. $n$
   1. Erzeugen Sie einen Vektor unabhängiger Standardnormalvariablen.$Z = (Z_1, \ldots, Z_d)'$
   2. Setze $X = AZ$
   3. Rückgabe .$U = (\Phi(X_1), \ldots, \Phi(X_d))'$

Der folgende Code in einer Beispielimplementierung dieses Algorithmus mit 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)))}))

Die folgende Tabelle zeigt die Daten, die sich aus dem obigen R-Code ergeben.