library(rgl)  #pro kresleni 3D grafu

# 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)  #nahodny vyber o rozsahu 1000 z dane kopuly 
plot(u)                        #graf realizace nahodneho vyberu
cor(u[,1],u[,2])               #vyberovy korelacni koeficient mezi slozkami
par(mfrow=c(1,2))
persp (indep.cop, dCopula)     #hustota dane kouply
contour(indep.cop, dCopula,nlevels=30)  #kontury hustoty dane kopuly


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

u <- rCopula(1000, Gumbel.cop)
plot(u)
cor(u[,1],u[,2]) 
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)
cor(u[,1],u[,2])
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)
cor(u[,1],u[,2])
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)
cor(u[,1],u[,2])
par(mfrow=c(1,2))
persp (norm.cop, dCopula)
contour(norm.cop, dCopula,nlevels=30)

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

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

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

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

u <- rCopula(1000, Galambos.cop)
plot(u)
cor(u[,1],u[,2])
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)
cor(u[,1],u[,2])
par(mfrow=c(1,2))
persp (Tawn.cop, dCopula)
contour(Tawn.cop, dCopula,nlevels=30)

##########################################
#Modelovani mnohorozmernych rozdeleni pomoci kopul
##########################################
#Chceme dvourozmerne rozdeleni, jehoz marginaly budou mit normalni a exponencialni
#rozdeleni. Zavislost budeme modelovat pomoci Gumbelovy kopuly.

normexp=mvdc(gumbelCopula(3), c("norm", "exp"),
list(list(mean= 0, sd = 1),list(rate = 1))) 

#hustota, distribucni funkce a kontury hustoty naseho rozdeleni:
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))


# Nahodny vyber z naseho dvourozmerneho rozdeleni
x.samp <- rMvdc(1000, normexp)
plot(x.samp)


# Odhady parametru modelu na zaklade nahodneho vyberu
fitMvdc(x.samp,normexp,c(1,1,1,5))

#pro srovnani uvedme odhady parametru jednotlivych slozek
x=x.samp[,1]
y=x.samp[,2]
mean(x)   # MLE pro mu
sd(x)     # MLE pro sigma
1/mean(y) # MLE pro lambda