library(rgl) 

# Soucinova kopula
x <- 1:1000/1000
y <- 1:1000/1000
z=matrix(rep(0,1000*1000),nrow=1000)

for (i in 1:1000){
for (j in 1:1000){
z[i,j]=x[i]*y[j]
}}

persp3d(x, y, z, col="skyblue",zlab='F(x,y)')
contour(x,y,z,xlab='x',ylab='y',main='Kontury F(x,y)')


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

persp3d(x, y, z, col="skyblue",zlab='F(x,y)')
contour(x,y,z,xlab='x',ylab='y',main='Kontury F(x,y)')


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

persp3d(x, y, z, col="skyblue",zlab='F(x,y)')
contour(x,y,z,xlab='x',ylab='y',main='Kontury F(x,y)')


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

persp3d(x, y, z, col="skyblue",zlab='F(x,y)')
contour(x,y,z,xlab='x',ylab='y',main='Kontury F(x,y)')


########################
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 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)')

contour(indep.cop, dCopula,nlevels=10,xlab='u',ylab='v',main='Kontury c(u,v)')
points(u,col=2,pch='+')

#(1) Gumbelova kopula

#(a)+(b)+(c)
theta=1.5
Gumbel.cop <- gumbelCopula(theta)

u <- rCopula(1000, Gumbel.cop)
plot(u,xlab='u',ylab='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)')

contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main='Kontury c(u,v)')
points(u,col=2,pch='+')

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

Gumbel.cop <- gumbelCopula(theta)
u <- rCopula(1000, Gumbel.cop)
plot(u,xlab='u',ylab='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)')

contour(Gumbel.cop, dCopula,nlevels=12,xlab='u',ylab='v',main='Kontury c(u,v)')
points(u,col=2,pch='+')

#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')))

#(3)

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

u <- rCopula(1000, Joe.cop)
plot(u,xlab='u',ylab='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)')


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

u <- rCopula(1000, Clayton.cop)
plot(u,xlab='u',ylab='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)')


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

u <- rCopula(1000, Frank.cop)
plot(u,xlab='x',ylab='y')
par(mfrow=c(1,2))
persp (Frank.cop, dCopula,xlab='u',ylab='v',main='f(u,v)')
contour(Frank.cop, dCopula,nlevels=30,xlab='x',ylab='y',main='Kontury c(u,v)')


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

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

u <- rCopula(1000, norm.cop)
plot(u,xlab='u',ylab='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)')

cor(u[,1],u[,2],method='kendall')
2/pi*asin(rho)
sin(pi/2*cor(u[,1],u[,2],method='kendall'))


#(2)
#Studentova t copula

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

u <- rCopula(1000, t.cop)
plot(u,xlab='u',ylab='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)')

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

#(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.5

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')


#(4)

################################################
#Vytvareni vicerozmernych rozdeleni pomoci kopul
################################################

#pomoci kopuly a marginalnich rozdeleni definujeme vlastni distribuci
mojerozdeleni=mvdc(gumbelCopula(5), c("unif", "gamma"),
list(list(min= 0, max = 100),list(shape=10, rate = 1/7)))

#nagenerujeme 1000 pozorovani z tohoto rozdeleni
u=rMvdc(1000, mojerozdeleni)
plot(u,xlab='x',ylab='y')

#vykreslime sdruzenou hustotu a distribucni funkci
persp  (mojerozdeleni, dMvdc, xlim = c(0, 100), ylim=c(30, 150),xlab='x', ylab='y', main = 'f(x,y)')
persp  (mojerozdeleni, pMvdc, xlim = c(0, 100), ylim=c(30, 150),xlab='x', ylab='y', main = 'F(x,y)')
contour(mojerozdeleni, dMvdc, xlim = c(0, 100), ylim=c(30, 150),nlevels=15,xlab='x', ylab='y',main='Kontury f(x,y)')

#pridejme jeste puvodni data
points(u,col=2,pch='+')