#clear
rm(list=ls())

##############################################################
#load packages / library
##############################################################

library(copula)
library(VineCopula)
# Fancy 3D plain scatterplots
library(scatterplot3d)
# ggplot2
library(ggplot2)
# Useful package to set ggplot plots one next to the other
library(grid)

##############################################################
#Import Data
##############################################################

mydata <- read.csv("C:/Users/Isariya/Desktop/mydata.csv")
x <- mydata[,2]
y <- mydata[,3]
plot(x,y)
plot(mydata$x,mydata$y)
hist(mydata$x)
hist(mydata$y)

##############################################################
#Parameter Estimation
##############################################################

#Estimate x  gamma distribution parameters and visually compare simulated vs observed data
x_mean <- mean(mydata$x)
x_var <- var(mydata$x)
x_rate <- x_mean / x_var
x_shape <- ( (x_mean)^2 ) / x_var

hist(mydata$x, breaks = 20, col = "green", density = 20)
hist(rgamma( nrow(mydata), rate = x_rate, shape = x_shape), breaks = 20,col = "blue", add = T, density = 20, angle = -45)

# Estimate y gamma distribution parameters and visually compare simulated vs observed data
y_mean <- mean(mydata$y)
y_var <- var(mydata$y)
y_rate <- y_mean / y_var
y_shape <- ( (y_mean)^2 ) / y_var

hist(mydata$y, breaks = 20, col = "green", density = 20)
hist(rgamma(nrow(mydata), rate = y_rate, shape = y_shape), breaks = 20, col = "blue", add = T, density = 20, angle = -45)

##############################################################
#correlation coefficients
##############################################################

?cor
cor(mydata, method = "pearson")
cor(mydata, method = "kendall")
cor(mydata, method = "spearman")

##############################################################
#fitting multivariate distribution with copula
##############################################################

?pobs
var_a <- pobs(mydata)[,2]
var_b <- pobs(mydata)[,3]
plot(var_a,var_b)

?BiCopSelect
selectedCopula <- BiCopSelect(var_a, var_b, familyset = NA)
selectedCopula

selectedCopula$family
selectedCopula$par

# Estimate copula parameters
cop_model <- claytonCopula(dim = 2)
m <- pobs(as.matrix(mydata[,2-3])
fit <- fitCopula(cop_model, m, method = 'ml')
coef(fit)