#########################################################
library(copula)
#####################################################
# Vypocet Kendallova tau a Spearmanova rho pro danou kopulu
gumbel.cop <- gumbelCopula(3)
tau(gumbel.cop)
rho(gumbel.cop)

# Vypocet parametru kopuly pro danou hodnotu Kendallova tau
tau=0.5
iTau(gumbelCopula(), tau)

# Umi pracovat i se 'slozitejsimi' kopulami
frank.cop <- frankCopula(3)
tau(frank.cop)
rho(frank.cop)

iTau(frankCopula(), tau)

########################
# Extreme value kopuly #
########################
# Galambosova copula
Galambos.cop <- galambosCopula(2.5)

u <- rCopula(1000, Galambos.cop)
plot(u,xlab='u',ylab='v')
par(mfrow=c(1,2))
persp (Galambos.cop, dCopula,xlab='u',ylab='v',main='c(u,v)')
contour(Galambos.cop, dCopula,nlevels=20,xlab='u',ylab='v',main='Kontury c(u,v)')

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


# Tawnova copula
Tawn.cop <- tawnCopula(c(0.93,0.1,5))

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

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

#####################################
# Nacteme soubor pojisteni.RData
#####################################
x=data[,1]  # Vek
y=data[,2]  # Vyse pojistneho plneni

plot(data,xlab='Vek',ylab='Vyse plneni',main='Data')

# Nejprve se podivame na obe promenne zvlast
hist(x,main='Histogram veku')
plot(density(x),main='Jadrovy odhad veku')

hist(y,main='Histogram vyse plneni')
plot(density(y),main='Jadrovy odhad vyse plneni')

# Definujeme dva modely (2 ruzne kopuly)

model1=mvdc(gumbelCopula(2), c("norm", "exp"),
list(list(mean= 40, sd = 10),list(rate = 1)))

model2=mvdc(normalCopula(0.5), c("norm", "exp"),
list(list(mean= 40, sd = 10),list(rate = 1/2)))


# Metoda maximalni verohodnosti
a=fitMvdc(data,model1,start=c(41,10,0.5,2))
summary(a)

# Ulozme si odhadnute parametry modelu
mu=coef(a)[[1]]
sd=coef(a)[[2]]
lambda=coef(a)[[3]]
theta=coef(a)[[4]]

# Porovnejme s odhady marginalnich rozdeleni
library(univariateML)
model_select(x)
model_select(y)

# Ulozme si je pro dalsi pouziti
mu_marginal=mean(x)
sd_marginal=sd(x)
lambda_marginal=1/mean(y)

# Vyzkousejme jeste druhy model

b=fitMvdc(data,model2,c(41,10,0.5,0.5))
summary(b)

# Zvolime prvni model

model_final=mvdc(gumbelCopula(theta), c("norm", "exp"),
list(list(mean= mu, sd = sd),list(rate = lambda)))

persp  (model_final, dMvdc, xlim = c(10, 70), ylim=c(0, 10),xlab='Vek', ylab='Vyse plneni', main = 'Sdruzena hustota')
persp  (model_final, pMvdc, xlim = c(10, 70), ylim=c(0, 10),xlab='Vek', ylab='Vyse plneni', main = 'Sdruzena distribucni funkce')
contour(model_final, dMvdc, xlim = c(10, 70), ylim=c(0, 10),nlevels=15,xlab='Vek', ylab='Vyse plneni',main='Kontury sdruzene hustoty')
points(data,col=2,pch='+')

# Metoda maximalni pseudo verohodnosti
u=pnorm(x,mu_marginal,sd_marginal)
v=pexp(y,lambda_marginal)

par(mfrow=c(1,2))
plot(data,xlab='Vek',ylab='Vyse plneni',main='Puvodni data')
plot(v~u,main='Transformovana data')

w=cbind(u,v) # Transformovana data

fitCopula(gumbelCopula(),w,method='mpl')

# Semiparametricky pristup
# Neparametricke odhady marginalnich rozdeleni
u2=pobs(data)
par(mfrow=c(1,2))
plot(data,xlab='Vek',ylab='Vyse plneni',main='Puvodni data')
plot(u2,main='Transformovana data',xlab='u',ylab='v')

fitCopula(gumbelCopula(),u2,method='mpl')

# Metoda inverze kandallova tau
fitCopula(gumbelCopula(),w,method='itau')
fitCopula(gumbelCopula(),u2,method='itau')

cor(x,y,method='kendall')
cor(u,v,method='kendall')
cor(u2[,1],u2[,2],method='kendall')

# Odhad parametru theta 'rucne'
1/(1-cor(x,y,method='kendall'))


# Vyber vhodne kopuly
library(VineCopula)
BiCopSelect(u2[,1], u2[,2])
BiCopSelect(u,v)


#####################
# Podminene hustoty
#####################
i=(0:1000)/100 # mrizka 0 az 10 po 0.01
par(mfrow=c(2,2))


x=30
fyx=function(y){
fyx=dMvdc(c(x,y),model_final)/dnorm(x,mu,sd)
return(fyx)
}
plot(Vectorize(fyx)(i)~i,type='l',main='Vek = 30',xlab='y',ylab='f(y|x=30)')

x=40
fyx=function(y){
fyx=dMvdc(c(x,y),model_final)/dnorm(x,mu,sd)
return(fyx)
}
plot(Vectorize(fyx)(i)~i,type='l',main='Vek = 40',xlab='y',ylab='f(y|x=40)')

x=50
fyx=function(y){
fyx=dMvdc(c(x,y),model_final)/dnorm(x,mu,sd)
return(fyx)
}
plot(Vectorize(fyx)(i)~i,type='l',main='Vek = 50',xlab='y',ylab='f(y|x=50)')

x=60
fyx=function(y){
fyx=dMvdc(c(x,y),model_final)/dnorm(x,mu,sd)
return(fyx)
}
plot(Vectorize(fyx)(i)~i,type='l',main='Vek = 60',xlab='y',ylab='f(y|x=60)')







 