# Example code to promote familiarity with the functionality of package EVcopula

####################################################################################################
# Preliminaries
####################################################################################################
# Installation
###############

####################################################################################################
# NOTE: The first set of instructions works fine at least on Ubuntu Linux; some variants may work on 
# Windows / Mac but see below for an alternative using the R library "devtools"
###################################################################################################

# Extract package from .tar.gz and set working directory to location where package is saved

setwd(...)

# Install using

install.packages("EVcopula",repos=NULL)

# To update to a future version, the above commands can be repeated. If the library is already loaded,
# then run these two commands first:

detach("package:EVcopula")
unloadNamespace("EVcopula")


# Alternative instructions:

library(devtools) # preceded by install.packages("devtools") if necessary
install_url("http://www.lancaster.ac.uk/~wadswojl/EVcopula.tar.gz")

####################################################################################################
# First loading
################

# After installation, load the library

library(EVcopula)

# For more info / to see a list of functions and their help files

?EVcopula


####################################################################################################
# Notes
################

# Use of the package "multicore" is recommended if you have >1 core, as computations can be slow. As
# a rough guide the first computation below took 43s using a single core and 33s using 4 cores.

# To use more than one core, use the argument "ncore" in key functions like fit.EV.copula to specify how 
# many you would like to use.

# Things will be slower the more uncensored datapoints there are, and typically when the dependence 
# is weaker.

####################################################################################################
# Example usage
################

# Bivariate Gaussian data, correlation rho

library(mvtnorm)

rho<-0.7
xy<-rmvnorm(1000,sigma=matrix(c(1,rho,rho,1),2,2))
uv<-cbind(unif(xy[,1]),unif(xy[,2]))

fit.Gauss<-fit.EV.copula(uv=uv,marg.thresh=0.95,norm="Lp",p=1)

# Same thing, but a little faster:
fit.Gauss<-fit.EV.copula(uv=uv,marg.thresh=0.95,norm="Lp",p=1,ncore=4)

# Different norm
fit.Gauss2<-fit.EV.copula(uv=uv,marg.thresh=0.95,norm="Linf",ncore=4)

# Check parameter estimates / convergence / AIC:

fit.Gauss
fit.Gauss2

# Diagnostics:

?plot.EV.copula

# just chi / chi bar (quick)

plot(fit.Gauss, which=c(1,2))

# check on the independence of the implied S, V (slower; see documentation and paper for definition of these)

plot(fit.Gauss, which=c(3))

# Inspect 2D likelihood surface:

# (Counter ranging from 0 to 1 is displayed to indicate progress)

Lik.Surface2D(fit.Gauss, lamrange = c(-0.25,0.4), arange = c(1,10), n=15)

# NB define this as an object if you wish to save the values for future use, e.g.

L2D.Gauss<-Lik.Surface2D(fit.Gauss, lamrange = c(-0.25,0.4), arange = c(1,10), n=15)

# Inspect profile likelihoods for the parameters (similar comments apply regarding the displayed counter and saving as an object)

Proflik1D(fit.Gauss,whichpar = "lambda",parrange = c(-0.25,0.41), npts = 15)
Proflik1D(fit.Gauss,whichpar = "alpha",parrange = c(1,10), npts = 15)

# Use fitted model to estimate a probability of interest:

Prob.Est(fit.Gauss, uvec=0.96,vvec=0.995)
Prob.Est(fit.Gauss, uvec=0.98,vvec=0.99)
Prob.Est(fit.Gauss, uvec=0.99,vvec=0.999)
Prob.Est(fit.Gauss, uvec=0.99999,vvec=0.9999)

Prob.Est(fit.Gauss, uvec=c(0.8,0.9999),vvec=c(0.995,0.99995))


# (Comparison of true probability in this case)

pmvnorm(lower=c(qnorm(0.96),qnorm(0.995)),sigma=matrix(c(1,rho,rho,1),2,2))
pmvnorm(lower=c(qnorm(0.98),qnorm(0.99)),sigma=matrix(c(1,rho,rho,1),2,2))
pmvnorm(lower=c(qnorm(0.99),qnorm(0.999)),sigma=matrix(c(1,rho,rho,1),2,2))
pmvnorm(lower=c(qnorm(0.99999),qnorm(0.9999)),sigma=matrix(c(1,rho,rho,1),2,2))

pmvnorm(lower=c(qnorm(0.8),qnorm(0.995)),sigma=matrix(c(1,rho,rho,1),2,2))-
  pmvnorm(lower=c(qnorm(0.8),qnorm(0.99995)),sigma=matrix(c(1,rho,rho,1),2,2))-
  pmvnorm(lower=c(qnorm(0.9999),qnorm(0.995)),sigma=matrix(c(1,rho,rho,1),2,2))+
  pmvnorm(lower=c(qnorm(0.9999),qnorm(0.99995)),sigma=matrix(c(1,rho,rho,1),2,2))

# Empirical

mean(uv[,1]>0.96&uv[,2]>0.995)
mean(uv[,1]>0.98&uv[,2]>0.99)
mean(uv[,1]>0.99&uv[,2]>0.999)
mean(uv[,1]>0.99999&uv[,2]>0.9999)

mean(uv[,1]>0.8&uv[,1]<0.9999&uv[,2]>0.995&uv[,2]<0.99995)

#################################################################################################

# Bivariate extreme value logistic data

library(evd)

xy<-rbvevd(1000,dep=0.5)
uv<-cbind(unif(xy[,1]),unif(xy[,2]))

fit.EVL<-fit.EV.copula(uv=uv,marg.thresh=0.95,norm="Lp",p=1,ncore=4)

# Different norm
fit.EVL2<-fit.EV.copula(uv=uv,marg.thresh=0.95,norm="Linf",ncore=4)

# Check parameter estimates / convergence / AIC:

fit.EVL
fit.EVL2

# Diagnostics:

?plot.EV.copula

# just chi / chi bar (quick)

plot(fit.EVL2, which=c(1,2))

# check on the independence of the implied S, V (slower; see documentation and paper for definition of these)

plot(fit.EVL2, which=c(3), ncore=4)

# Inspect 2D likelihood surface:

Lik.Surface2D(fit.EVL2, lamrange = c(0.2,2), arange = c(0.1,5), n=15, ncore=4)

# Profile likelihood for lambda

Proflik1D(fit.EVL2,whichpar = "lambda",parrange = c(0.5,2), npts = 15)

# Use fitted model to estimate a probability of interest:

Prob.Est(fit.EVL2, uvec=0.96,vvec=0.995)
Prob.Est(fit.EVL2, uvec=0.98,vvec=0.99)
Prob.Est(fit.EVL2, uvec=0.99,vvec=0.999)
Prob.Est(fit.EVL2, uvec=0.99999,vvec=0.9999)

Prob.Est(fit.EVL2, uvec=c(0.8,0.9999),vvec=c(0.995,0.99995))


# (Comparison of true probability in this case)

qgumb<-function(u){-log(-log(u))}

pbvevd(q=c(qgumb(0.96),qgumb(0.995)),dep=0.5,low=F)
pbvevd(q=c(qgumb(0.98),qgumb(0.99)),dep=0.5,low=F)
pbvevd(q=c(qgumb(0.99),qgumb(0.999)),dep=0.5,low=F)
pbvevd(q=c(qgumb(0.99999),qgumb(0.9999)),dep=0.5,low=F)

pbvevd(q=c(qgumb(0.8),qgumb(0.995)),dep=0.5,low=F)-
  pbvevd(q=c(qgumb(0.8),qgumb(0.99995)),dep=0.5,low=F)-
  pbvevd(q=c(qgumb(0.9999),qgumb(0.995)),dep=0.5,low=F)+
  pbvevd(q=c(qgumb(0.9999),qgumb(0.99995)),dep=0.5,low=F)

# Empirical

mean(uv[,1]>0.96&uv[,2]>0.995)
mean(uv[,1]>0.98&uv[,2]>0.99)
mean(uv[,1]>0.99&uv[,2]>0.999)
mean(uv[,1]>0.99999&uv[,2]>0.9999)

mean(uv[,1]>0.8&uv[,1]<0.9999&uv[,2]>0.995&uv[,2]<0.99995)