mick001

# Estimate copula parameters
cop_model <- claytonCopula(dim = 2)
m <- pobs(as.matrix(mydata))
fit <- fitCopula(cop_model, m, method = 'ml')
coef(fit)
##    param 
## 1.482168 

# Check Kendall's tau value for the Clayton copula with theta = 1.48
tau(claytonCopula(param = 1.48))

mick001

var_a <- pobs(mydata)[,1]
var_b <- pobs(mydata)[,2]
selectedCopula <- BiCopSelect(var_a, var_b, familyset = NA)
selectedCopula

## $p.value.indeptest
## [1] NA
## 
## $family
## [1] 3

mick001

# Measure association using Kendall's Tau
cor(mydata, method = "kendall")

##           x         y
## x 1.0000000 0.4212052
## y 0.4212052 1.0000000

mick001

# Estimate x  gamma distribution parameters and visually compare simulated vs observed data
x_mean <- mean(mydata$x)
x_var <- var(mydata$x)
x_rate <- x_mean / x_var
x_shape <- ( (x_mean)^2 ) / x_var

hist(mydata$x, breaks = 20, col = "green", density = 20)
hist(rgamma( nrow(mydata), rate = x_rate, shape = x_shape), breaks = 20,col = "blue", add = T, density = 20, angle = -45)

# Estimate y gamma distribution parameters and visually compare simulated vs observed data

mick001

# Copula package
library(copula)
# Fancy 3D plain scatterplots
library(scatterplot3d)
# ggplot2
library(ggplot2)
# Useful package to set ggplot plots one next to the other
library(grid)
set.seed(235)

mick001

# Generate a bivariate normal copula with rho = 0.7
normal <- normalCopula(param = 0.7, dim = 2)
# Generate a bivariate t-copula with rho = 0.8 and df = 2
stc <- tCopula(param = 0.8, dim = 2, df = 2)

mick001

# Build a Frank, a Gumbel and a Clayton copula
frank <- frankCopula(dim = 2, param = 8)
gumbel <- gumbelCopula(dim = 3, param = 5.6)
clayton <- claytonCopula(dim = 4, param = 19)

# Print information on the Frank copula
print(frank)

mick001

# Select the copula
cp <- claytonCopula(param = c(3.4), dim = 2)

# Generate the multivariate distribution (in this case it is just bivariate) with normal and t marginals
multivariate_dist <- mvdc(copula = cp,
              margins = c("norm", "t"),
              paramMargins = list(list(mean = 2, sd=3),
                                  list(df = 2)) )

print(multivariate_dist)

mick001

# Generate random samples
fr <- rCopula(2000, frank)
gu <- rCopula(2000, gumbel)
cl <- rCopula(2000, clayton)

# Plot the samples
p1 <- qplot(fr[,1], fr[,2], colour = fr[,1], main="Frank copula random samples theta = 8", xlab = "u", ylab = "v")
p2 <- qplot(gu[,1], gu[,2], colour = gu[,1], main="Gumbel copula random samples theta = 5.6", xlab = "u", ylab = "v") 
p3 <- qplot(cl[,1], cl[,2], colour = cl[,1], main="Clayton copula random samples theta = 19", xlab = "u", ylab = "v")

mick001

# Generate the normal copula and sample some observations
coef_ <- 0.8
mycopula <- normalCopula(coef_, dim = 2)
u <- rCopula(2000, mycopula)

# Compute the density
pdf_ <- dCopula(u, mycopula)

# Compute the CDF
cdf <- pCopula(u, mycopula)