Simulating from a Bivariate Gaussian Copula

The 2-dimensional Gaussian Copula is a distribution on the square \([0,1]\times [0,1]\). For a given correlation matrix \(R = \begin{pmatrix} 1 & \rho \\ \rho & 1\end{pmatrix}\), the Gaussian copula cumulative distribution function (cdf) \[C_{R}(u_1,u_2)=\Phi _{R}\left(\Phi ^{-1}(u_{1}),\Phi ^{-1}(u_{2})\right),\] where \(\Phi_R\) is the bivariate normal cdf with covariance matrix \(R\) (equal to the correlation matrix), and \(\Phi^{-1}\) is the quantile normal function (the inverse of the standard normal cdf). Simulating from this distribution can be done as follows

Simulate \((x_i,y_i)^{\top} \sim N_2(0,R)\), \(i=1,\dots,n\).
Calculate \((v_i,w_i)^{\top} = \left(\Phi(x_i),\Phi(y_i)\right)^{\top}\).

Then, the sample \((v_i,w_i)^{\top}\), \(i=1,\dots,n\), is a sample from a bivariate Gaussian copula with association parameter \(\rho\).

Simulating from a bivariate distribution, based on the Gaussian copula, with given marginals

Suppose now that we are interested on simulating from a bivariate distribution with marginals \(F_1\) and \(F_2\). A possible way for doing so consists of constructing a bivariate distribution using a Gaussian copula. Simulating from such distribution can be done as follows

Simulate \((x_i,y_i)^{\top} \sim N_2(0,R)\), \(i=1,\dots,n\).
Calculate \((v_i,w_i)^{\top} = \left(F_1^{-1}(\Phi(x_i)),F_2^{-1}(\Phi(y_i))\right)^{\top}\).

The following R code shows how to simulate from a bivariate distribution with marginals \(F_1\) and \(F_2\), based on a Gaussian copula.

Other things to consider

There exist many (infinite) types of copulas.
Other copula functions can be used to capture different types of dependencies. One example of this is the t-copula.
You can fit a bivariate distribution that was obtained by using a Gaussian (or a t-copula) using the maximum likelihood estimation method.
The ideas presented here can be easily extended to higher dimensions, but we have focused on the bivariate case for illustrative purposes.

R code

rm(list=ls())
# Required packages
library(mvtnorm)
library(ggplot2)
library(ggExtra)

## Warning: package 'ggExtra' was built under R version 3.5.3

library(sn)

## Loading required package: stats4

## 
## Attaching package: 'sn'

## The following object is masked from 'package:stats':
## 
##     sd

library(twopiece)

## Loading required package: flexclust

## Loading required package: grid

## Loading required package: lattice

## Loading required package: modeltools

## Loading required package: label.switching

##################################################################################################
# Function to simulate from a bivariate distribution with (possible) different marginals 
# using a bivariate Gaussian copula
##################################################################################################
# n: sample size
# qmarg1: quantile function for marginal 1
# qmarg2: quantile function for marginal 2
# rho: correlation parameter in Gaussian Copula


sim.GC <- function(n, rho, qmarg1, qmarg2){
R <- rbind(c(1,rho),c(rho,1))
dat <- rmvnorm(n, mean = c(0,0), sigma = R)
dat[,1] <- qmarg1(pnorm(dat[,1]))
dat[,2] <- qmarg2(pnorm(dat[,2]))
return(dat)
}

##################################################################################################
# Examples
##################################################################################################
set.seed(1234)
#--------------------
# Uniform marginals
#--------------------
q1 <- function(p) qunif(p)
q2 <- function(p) qunif(p)

ns <- 1000; rho = 0.5
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.5")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.75
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.75")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.95
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.95")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()


#--------------------
# Logistic marginals
#--------------------
q1 <- function(p) qlogis(p)
q2 <- function(p) qlogis(p)

ns <- 1000; rho = 0.5
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.5")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.75
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.75")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.95
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.95")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

#----------------------------------------------------
# Student-t with 4 degrees of freedom marginals
#----------------------------------------------------
q1 <- function(p) qt(p, df = 4)
q2 <- function(p) qt(p, df = 4)

ns <- 1000; rho = 0.5
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.5")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.75
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.75")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.95
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.95")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()


#----------------------------------------------------
# Student-t with 4 degrees of freedom marginal for X
# Standard normal for Y
#----------------------------------------------------
q1 <- function(p) qt(p, df = 4)
q2 <- function(p) qnorm(p)

ns <- 1000; rho = 0.5
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.5")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.75
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.75")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.95
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.95")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()


#----------------------------------------------------
# Skew-normal marginals
#----------------------------------------------------
q1 <- function(p) qsn(p,0,1,5)
q2 <- function(p) qsn(p,0,1,-5)

ns <- 1000; rho = 0.5
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.5")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.75
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.75")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.95
sim <- sim.GC(n = ns,rho = rho,q1,q2)
df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.95")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()



#----------------------------------------------------
# twopiece-normal marginals
#----------------------------------------------------
q1 <- function(p) qtp3(p,0,1,0.5, param= "eps", FUN = qnorm)
q2 <- function(p) qtp3(p,0,1,-0.5, param= "eps", FUN = qnorm)

ns <- 1000; rho = 0.5
sim <- sim.GC(n = ns,rho = rho,q1,q2)

## Warning in FUN((p - gamma)/(1 - gamma)): NaNs produced

## Warning in FUN(p/(1 + gamma)): NaNs produced

df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.5")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.75
sim <- sim.GC(n = ns,rho = rho,q1,q2)

## Warning in FUN((p - gamma)/(1 - gamma)): NaNs produced

## Warning in FUN((p - gamma)/(1 - gamma)): NaNs produced

df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.75")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

ns <- 1000; rho = 0.95
sim <- sim.GC(n = ns,rho = rho,q1,q2)

## Warning in FUN((p - gamma)/(1 - gamma)): NaNs produced

## Warning in FUN((p - gamma)/(1 - gamma)): NaNs produced

df <- data.frame(x = sim[,1], y = sim[,2])
p <- ggplot(df, aes(x, y)) + geom_point() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) +
  xlab("X") + ylab("Y") + ggtitle(expression(paste(rho, "=0.95")))
ggExtra::ggMarginal(p, type = "histogram")

grid::grid.newpage()

Simulating from a Bivariate Gaussian Copula

F. Javier Rubio

October 21, 2019

Simulating from a Bivariate Gaussian Copula

Simulating from a bivariate distribution, based on the Gaussian copula, with given marginals

Other things to consider

R code