library(VineCopula)
library(copula)

u <- pobs(as.matrix(cbind(stock1,stock2)))[,1]
v <- pobs(as.matrix(cbind(stock1,stock2)))[,2]
selectedCopula <- BiCopSelect(u,v,familyset=NA)

Bivariate copula: BB8 (par = 1.78, par2 = 0.86, tau = 0.19)

copula <- BB8Copula(param = c(1.78,0.86))
m <- pobs(as.matrix(cbind(stock1,stock2)))
fit <- fitCopula(copula,m,method='itau')
coef(fit)

my_dist <- mvdc(claytonCopula(param = 1.48, dim = 2), margins = c("gamma","gamma"), paramMargins = list(list(shape = x_shape, rate = x_rate), list(shape = y_shape, rate = y_rate)))

sim <- rMvdc(306, my_dist)
cor(cbind(stock1,stock2), method = "kendall")
cor(sim, method = "kendall")

> cor(cbind(stock1,stock2), method = "kendall")
        stock1     stock2
stock1 1.0000000 0.1955256
stock2 0.1955256 1.0000000
> cor(sim, method = "kendall")
          [,1]      [,2]
[1,] 1.0000000 0.4441659
[2,] 0.4441659 1.0000000