library(rgl) 

# Soucinova kopula
u <- 1:1000/1000
v <- 1:1000/1000
C=matrix(rep(0,1000*1000),nrow=1000)

for (i in 1:1000){
for (j in 1:1000){
C[i,j]<-u[i]*v[j]
}}

persp3d(u, v, C, col="skyblue",zlab='C(u,v)')
contour(u,v,C,xlab='u',ylab='v',main='Kontury C(u,v)')


# Horni kopula
for (i in 1:1000){
for (j in 1:1000){
C[i,j]<-min(u[i],v[j])
}}

persp3d(u, v, C, col="skyblue",zlab='C(u,v)')
contour(u,v,C,xlab='u',ylab='v',main='Kontury C(u,v)')


# Dolni kopula
for (i in 1:1000){
for (j in 1:1000){
C[i,j]<-max(u[i]+v[j]-1,0)
}}

persp3d(u, v, C, col="skyblue",zlab='C(u,v)')
contour(u,v,C,xlab='u',ylab='v',main='Kontury C(u,v)')


# Gumbelova kopula
theta<-2
for (i in 1:1000){
for (j in 1:1000){
C[i,j]=exp(-((-log(u[i]))^(theta)+(-log(v[j]))^(theta))^(1/theta) )
}}

persp3d(u, v, C, col="skyblue",zlab='C(u,v)')
contour(u,v,C,xlab='u',ylab='v',main='Kontury C(u,v)')



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

# Kopula nezavislosti
indep.cop <- indepCopula(2)

#nagenerujeme 1000 pozorovani z teto kopuly
u<-rCopula(1000, indep.cop)
plot(u,xlab='u',ylab='v')

#vykreslime kopulu a jeji kontury
par(mfrow=c(1,2))
persp (indep.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
contour(indep.cop, pCopula,nlevels=10,xlab='u',ylab='v',main='Kontury C(u,v)')

#vykreslime jeste hustotu kopuly a jeji kontury
par(mfrow=c(1,2))
persp (indep.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(indep.cop, dCopula,nlevels=10,xlab='u',ylab='v',main='Kontury c(u,v)')


#(1) Gumbelova kopula

#(a)
theta<-1.5
Gumbel.cop <- gumbelCopula(theta)

persp (Gumbel.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (Gumbel.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main='Kontury c(u,v)')

#(b)
u <- rCopula(1000, Gumbel.cop)
#(c)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main='Kontury c(u,v)',add=TRUE)

#(d)
# pro danou hodnotu Kendallova tau urcime parametr koupuly theta
tau<-0.5
theta<-1/(1-tau)
theta

Gumbel.cop <- gumbelCopula(theta)

persp (Gumbel.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (Gumbel.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main='Kontury c(u,v)')

u <- rCopula(1000, Gumbel.cop)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main='Kontury c(u,v)',add=TRUE)

#overme jeste empirickym odhadem
cor(u[,1],u[,2],method='kendall')

#(e)
#zavislost tvaru kopuly na parametru theta

par(mfrow=c(2,2))
theta=1.1
Gumbel.cop <- gumbelCopula(theta)
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main=expression(paste(theta,'=1.1')))

theta=2
Gumbel.cop <- gumbelCopula(theta)
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main=expression(paste(theta,'=2')))

theta=5
Gumbel.cop <- gumbelCopula(theta)
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main=expression(paste(theta,'=5')))

theta=20
Gumbel.cop <- gumbelCopula(theta)
contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main=expression(paste(theta,'=20')))

#(2)
#(a)
rho<--0.7
nu<-5
t.cop <- tCopula(rho,df=nu)

persp (t.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (t.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main='Kontury c(u,v)')

#(b)
u <- rCopula(1000, t.cop)

#(c)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main='Kontury c(u,v)',add=TRUE)

#(d)
#zavislost tvaru kopuly na parametru rho

par(mfrow=c(2,2))
nu<-5
rho<--0.8

t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(rho,'=-0.8')))

rho<--0.1
t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(rho,'=-0.1')))

rho<-0.1
t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(rho,'=0.1')))

rho<-0.8
t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(rho,'=0.8')))

#(e)
#zavislost tvaru kopuly na parametru nu

par(mfrow=c(2,2))
nu<-1
rho<--0.7

t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(nu,'=1')))

nu<-2
t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(nu,'=2')))

nu<-5
t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(nu,'=5')))

nu<-20
t.cop <- tCopula(rho,df=nu)
contour(t.cop, dCopula,nlevels=30,xlab='u',ylab='v',main=expression(paste(nu,'=20')))

#(f)
#teoreticka hodnota Kendallova tau
2/pi*asin(rho)

#empiricka hodnota Kendallova tau
cor(u[,1],u[,2],method='kendall')



#(3)

# Joeova kopula
theta<-1.5
Joe.cop <- joeCopula(theta)

persp (Joe.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (Joe.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(Joe.cop, dCopula,nlevels=15,xlab='u',ylab='v',main='Kontury c(u,v)')

u <- rCopula(1000, Joe.cop)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(Joe.cop, dCopula,nlevels=15,xlab='u',ylab='v',main='Kontury c(u,v)',add=TRUE)


# Claytonova kopula
theta<-5
Clayton.cop <- claytonCopula(theta)

persp (Clayton.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (Clayton.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(Clayton.cop, dCopula,nlevels=30,xlab='u',ylab='v',main='Kontury c(u,v)')

u <- rCopula(1000, Clayton.cop)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(Clayton.cop, dCopula,nlevels=30,xlab='u',ylab='v',main='Kontury c(u,v)',add=TRUE)


# Frankova kopula
theta<-4
Frank.cop <- frankCopula(theta)

persp (Frank.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (Frank.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(Frank.cop, dCopula,nlevels=30,xlab='x',ylab='y',main='Kontury c(u,v)')

u <- rCopula(1000, Frank.cop)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(Frank.cop, dCopula,nlevels=30,xlab='x',ylab='y',main='Kontury c(u,v)',add=TRUE)


# Normalni kopula
rho<-0.6
norm.cop <- normalCopula(rho)

persp (norm.cop, pCopula,xlab='u',ylab='v',main='C(u,v)')
par(mfrow=c(1,2))
persp (norm.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(norm.cop, dCopula,nlevels=30,xlab='u',ylab='v',main='Kontury c(u,v)')

u <- rCopula(1000, norm.cop)
plot(u,xlab='u',ylab='v',col=2,pch='+')
contour(norm.cop, dCopula,nlevels=30,xlab='u',ylab='v',main='Kontury c(u,v)',add=TRUE)


#(4)
#pomoci kopuly a marginalnich rozdeleni definujeme vlastni distribuci

mojerozdeleni<-mvdc(gumbelCopula(3), c("unif", "gamma"),
list(list(min=0, max=100),list(shape=10, rate=1/7)))

#(a)
persp  (mojerozdeleni, pMvdc, xlim = c(0, 100), ylim=c(30, 150),xlab='x', ylab='y', main = 'F(x,y)')

#(b)
persp  (mojerozdeleni, dMvdc, xlim = c(0, 100), ylim=c(30, 150),xlab='x', ylab='y', main = 'f(x,y)')

#(c)
contour(mojerozdeleni, dMvdc, xlim = c(0, 100), ylim=c(30, 150),nlevels=15,xlab='x', ylab='y',main='Kontury f(x,y)')

#(d)
u<-rMvdc(1000, mojerozdeleni)
points(u,col=2,pch='+')