library(rgl) 

# Soucinova kopula
x <- 1:1000/1000
y <- 1:1000/1000
z <- x %o% y

persp3d(x, y, z, col="skyblue")
contour(x,y,z)


# Horni kopula
for (i in 1:1000){
for (j in 1:1000){
u=x[i]
v=y[j]
z[i,j]=min(u,v)
}}

persp3d(x, y, z, col="skyblue")
contour(x,y,z)


# Dolni kopula
for (i in 1:1000){
for (j in 1:1000){
u=x[i]
v=y[j]
z[i,j]=max(u+v-1,0)
}}

persp3d(x, y, z, col="skyblue")
contour(x,y,z)


# Gumbelova kopula
theta=3
for (i in 1:1000){
for (j in 1:1000){
u=x[i]
v=y[j]
z[i,j]=exp(-((-log(u))^(theta)+(-log(v))^(theta))^(1/theta) )
}}

persp3d(x, y, z, col="skyblue")
contour(x,y,z)


#########################################################
library(copula)
########################
# Archimedovske kopuly #
########################

# Independence copula
indep.cop <- indepCopula(2)

u <- rCopula(1000, indep.cop)
plot(u)
par(mfrow=c(1,2))
persp (indep.cop, dCopula)
contour(indep.cop, dCopula,nlevels=30)


# Gumbel copula
Gumbel.cop <- gumbelCopula(3)

u <- rCopula(1000, Gumbel.cop)
plot(u)
par(mfrow=c(1,2))
persp (Gumbel.cop, dCopula)
contour(Gumbel.cop, dCopula,nlevels=30)


# Joe copula
Joe.cop <- joeCopula(2)

u <- rCopula(1000, Joe.cop)
plot(u)
par(mfrow=c(1,2))
persp (Joe.cop, dCopula)
contour(Joe.cop, dCopula,nlevels=30)


# AMH copula
AMH.cop <- amhCopula(0.9)

u <- rCopula(1000, AMH.cop)
plot(u)
par(mfrow=c(1,2))
persp (AMH.cop, dCopula)
contour(AMH.cop, dCopula,nlevels=30)

####################
# Elipticke kopuly #
####################

# Gausian copula
norm.cop <- normalCopula(0.6)

u <- rCopula(1000, norm.cop)
plot(u)
par(mfrow=c(1,2))
persp (norm.cop, dCopula)
contour(norm.cop, dCopula,nlevels=30)

# t copula
t.cop <- tCopula(0.1,df=5)

u <- rCopula(1000, t.cop)
plot(u)
par(mfrow=c(1,2))
persp (t.cop, dCopula)
contour(t.cop, dCopula,nlevels=30)

########################
# Extreme value kopuly #
########################

# Galambos copula
Galambos.cop <- galambosCopula(2.5)

u <- rCopula(1000, Galambos.cop)
plot(u)
par(mfrow=c(1,2))
persp (Galambos.cop, dCopula)
contour(Galambos.cop, dCopula,nlevels=30)

# Tawn copula
Tawn.cop <- tawnCopula(c(0.93,0.5,0.2))

u <- rCopula(1000, Tawn.cop)
plot(u)
par(mfrow=c(1,2))
persp (Tawn.cop, dCopula)
contour(Tawn.cop, dCopula,nlevels=30)

##########################################
#Multivariate distributions via copulas
##########################################
normexp=mvdc(gumbelCopula(3), c("norm", "exp"),
list(list(mean= 0, sd = 1),list(rate = 1)))

persp  (normexp, dMvdc, xlim = c(-3,3), ylim=c(0, 8), main = "N(0,1) + Exp(1)")
persp  (normexp, pMvdc, xlim = c(-3, 3), ylim=c(0, 8), main = "N(0,1) + Exp(1)")
contour(normexp, dMvdc, xlim = c(-3, 3), ylim=c(0, 8))


# Generate (bivariate) random numbers from that, and visualize
x.samp <- rMvdc(1000, normexp)
plot(x.samp)


# Fit a model
fitMvdc(x.samp,normexp,c(1,1,1,5))
x=x.samp[,1]
y=x.samp[,2]
mean(x)   # MLE pro mu
sd(x)     # MLE pro sigma
1/mean(y) # MLE pro lambda