library(copula)
library(MASS)
library(kde)
library(sm)


# 1. Generate synthetic data, standard normal marginals
# Read data - Lees Ferry monthly flows

flows=matrix(scan("http://civil.colorado.edu/~balajir/CVEN5454/r-project-data/LeesFerry-monflow-1906-2006.txt"),ncol=13,byrow=T)
flows=flows[,2:13]/10^6

## Select Apr and May flows. May flows has bimodality
data = flows[,4:5]
u = pobs(data)

## fit Gamma for Apr flows
zgamma = fitdistr(data[,1], "gamma")
uu = pgamma(data[,1], shape=zgamma$estimate[1], scale=1/zgamma$estimate[2])

## Gamma CDF for April and Empirical CDF for May
u = cbind(uu,pobs(data)[,2])


# 2. Fit the copula 
fnc = fitCopula(normalCopula(dim=2,disp='un'),u)

# 3. random generate from the fitted copula
cop_sim = rCopula(10000,normalCopula(fnc@estimate,dim=2,dispstr='un'))

# 4.  convert back to data space, assuming standard normal marginals

## X comes from Gamma and Y from the sample
data_sim_x = qgamma(cop_sim[,1],shape=zgamma$estimate[1], scale=1/zgamma$estimate[2])
data_sim_y = quantile(data[,2],cop_sim[,2])

# 5. plot contours of the simulated data and the data together
contour(kde2d(data_sim_x, data_sim_y,n=100))
points(data[,1],data[,2],cex=0.5,pch=19)


contour(kde2d(data_sim_x, data_sim_y,n=100),nlevels=20)

or specify the contour levels:

contour(kde2d(data_sim_x, data_sim_y,n=100),levels=c(0.3,0.2,0.1,0.05,0.01,0.001))

## Plot the May histogram

hist(data[,2],probability=TRUE)
N = dim(data)[1]

fc = frankCopula(dim=2)
ffc=fitCopula(fc,u)

### Simulate
for(isim in 1:100){


zsim = rCopula(N,frankCopula(coef(ffc),dim=2))

## Get the Kernel quantile..
fhat = kde(data[,2])
simmay=qkde(zsim[,2],fhat)

sm.density(simmay,add=TRUE,col="grey")

}

sm.density(data[,2],add=TRUE,col="red",lwd=2)