# R/Densities.R In ExtremalDep: Extremal Dependence Models

#### Documented in angulardens

```################################################################################################
### Authors: Boris Beranger and Simone Padoan        	 		                    							 ###
### 																							                                           ###
###	Emails: [email protected], [email protected]				            				 ###
### 																						                                          	 ###
###	Institutions: Department of Decision Sciences, University Bocconi of Milan		      		 ###
### School of Mathematics and Statistics, University of New South Wales 			        			 ###
### 																							                                           ###
### File name: Densities.r	                 							                            	     ###
### 																						                                          	 ###
### Description:                                  		              					      		     ###
### This file provides the Angular densities of extremal dependence models: 					       ###
### 1) Pairwise Beta, Dirichlet and Husler-Reiss (mass only on the interior of the simplex)  ###
### 2) Asymmetric logistic and Extremal-t (mass on all subsets of the simplex)				       ###
### It also provides the likelihood and log-likelihood functions								             ###
### 																							                                           ###
### Last change: 06/08/2019                         		  									                 ###
### 																					                                          		 ###
################################################################################################

dens <- function(x=rbind(c(0.1,0.3,0.6),c(0.1,0.2,0.7)), model="Pairwise", par=c(2,2,2,4), c=NULL, log=FALSE, vectorial=TRUE){

### Functions needed for the Pairwise Beta model

# Model defined for a minimum of three dimensions

dens_pb <- function (x, b, alpha, log, vectorial){
if(any(b<=0) || b<=0){return(1e-300)}
hij<-function(wi,wj,bij,alpha,p){
wij=wi+wj;
return(wij^(2*alpha-1)*(1-wij)^(alpha*(p-2)-p+2)*gamma(2*bij)/gamma(bij)^2*(wi/wij)^(bij-1)*(wj/wij)^(bij-1))
}

dens_pairb<-function(w,b,alpha){
p<-length(w);
dens=0;
K= 2*factorial(p-3)*gamma(alpha*p+1)/(p*(p-1)*gamma(2*alpha+1)*gamma(alpha*(p-2)))
for(i in 1:(p-1)){
for(j in (i+1):p){
d= hij(w[i],w[j],b[(i-1)*p-sum(1:i)+j],alpha,p)
dens=dens+d
}
}
dens=dens*K
return(dens)
}

xvect = as.double(as.vector(t(x)))
if (is.vector(x)) {
dim = as.integer(length(x))
n = as.integer(1)
if(round(sum(x),7) != 1){ stop("Data is not angular")}
result <- dens_pairb(x,b,alpha)
}
else {
dim = as.integer(ncol(x))
n = as.integer(nrow(x))
if (sum(apply(x,1,sum)) != n){ stop("Data is not angular")}

if (vectorial) {
result = double(n)
result <- apply(x,1,function(y){dens_pairb(y,b,alpha)})
} else { # vectorial=FALSE mean we return the likelihood function
result = as.double(1)
result <- prod(apply(x,1,function(y){dens_pairb(y,b,alpha)}))
}
}
if(log)
return(log(result))
else return(result)
}

### Functions needed for the Husler-Reiss model

dens_hr <- function (x, lambda, log, vectorial){

dens_husler<-function(w,lambda){
p<-length(w);
k=p-1
w.tilde<-rep(0,k);
Sigma<-diag(k);

for(i in 1:k){
if(w[1]==0 | w[i+1]==0){return(1e-50)} else{
w.tilde[i] <- log(w[i+1]/w[1])+2*lambda[i]^2;

for(j in i:k){
if(i==j){Sigma[i,j]=Sigma[j,i]=2*(lambda[i]^2+lambda[j]^2)}else{
Sigma[i,j]=Sigma[j,i]=2*(lambda[i]^2+lambda[j]^2-lambda[sum(k:(k-i+1))+j-i]^2)
}
}
}
}

if(sum(eigen(Sigma)\$values<0)>=1){return(1e-50)} # Check if matrix is positive definite
if(any(is.na(w.tilde))){return (1e-50)}

part1 = w[1]^2*prod(w[2:p])*(2*pi)^(k/2)* abs(det(Sigma))^(1/2)
part2 = exp(-0.5*t(w.tilde) %*% solve(Sigma) %*% t(t(w.tilde)))

return(part2/part1)
}

xvect = as.double(as.vector(t(x)))
if (is.vector(x)) {
dim = as.integer(length(x))
n = as.integer(1)
if(round(sum(x),7) != 1){ stop("Data is not angular")}
if(log){
result = log(dens_husler(x,lambda))
}else{
result = dens_husler(x,lambda)
}
}
else {
dim = as.integer(ncol(x))
n = as.integer(nrow(x))
if (sum(apply(x,1,sum)) != n){ stop("Data is not angular")}

if (vectorial) {
result = double(n)
if (log){
result <- apply(x,1,function(y){log(dens_husler(y,lambda))})
}else{
result <- apply(x,1,function(y){dens_husler(y,lambda)})
}
} else { # vectorial=FALSE mean we return the likelihood function
result = as.double(1)
if (log){
result <- sum(apply(x,1,function(y){log(dens_husler(y,lambda))}))
}else{
result <- prod(apply(x,1,function(y){dens_husler(y,lambda)}))
}
}
}
return(result)
}

### Functions needed for the Dirichlet model

dens_di <- function (x, para, log, vectorial){

dens_diri <- function(w,para){
d <- length(para)

if(any(para <=0)){return(1e-50)}
if(length(para) != d){stop("Wrong length of parameter")}

part1 <- prod(para/gamma(para))
part2 <- gamma(sum(para)+1)/(sum(para*w)^(d+1))
part3 <- prod((para*w/sum(para*w))^(para-1))
return(part1*part2*part3)
}

xvect = as.double(as.vector(t(x)))
if (is.vector(x)) {
dim = as.integer(length(x))
n = as.integer(1)
if(round(sum(x),7) != 1){ stop("Data is not angular")}
result = dens_diri(x,para)
}
else {
dim = as.integer(ncol(x))
n = as.integer(nrow(x))
if (sum(apply(x,1,sum)) != n){ stop("Data is not angular")}

if (vectorial) {
result = double(n)
result <- apply(x,1,function(y){dens_diri(y,para)})
} else { # vectorial=FALSE mean we return the likelihood function
result = as.double(1)
result <- prod(apply(x,1,function(y){dens_diri(y,para)}))
}
}
if(log)
return(log(result))
else return(result)
}

### Functions needed for the Asymmetric logistic model

dens_al <- function (x, alpha, beta, c, log, vectorial){

# The 2d case:

# in the following functions:
#
# alpha is a vector of size 1: for the subset {1,2}
# beta is a vector of size 2: for [1,{1,2}] and [2,{1,2}]
# beta for [1,{1}], [2,{2}] and [3,{3}] are omitted as obtained as [1,{1}] = 1 - [1,{1,2}] and [2,{2}] = 1 - [2,{1,2}]

interior_alm_2d <- function(w,alpha,beta){
part1 <- (alpha-1) * (beta[1]*beta[2])^alpha * (w[1]*w[2])^(alpha-2)
part2 <- ( (beta[1]*w[2])^alpha + (beta[2]*w[2])^alpha )^(1/alpha-2)

return(part1*part2)
}

corners_alm_2d <- function(w,alpha,beta,s){ # mass on the corner of the s-th component
if(s==1){return(1-beta[1])}
if(s==2){return(1-beta[2])}
}

dens_alm_2d <- function(w,alpha,beta,c){
if(length(alpha)!=1){return(stop("Wrong length of parameter alpha"))}
if(length(beta)!=2){return(stop("Wrong length of parameter beta"))}
if( (alpha < 1) || any(beta < 0) || any(beta > 1)){return(1e-50)}

if(sum(w > (1-c) ) == 1){ # then we are in a corner
ind <- which(w > (1-c))
return(corners_alm_2d(w,alpha,beta,ind))
} else {
return(interior_alm_2d(w,alpha,beta))
}
}

# The 3d case:

# in the following functions:
#
# alpha is a vector of size 4: for the subsets {1,2}, {1,3}, {2,3} and {1,2,3}
# beta is a vector of size 9: for [1,{1,2}], [2,{1,2}], [1,{1,3}], [3,{1,3}], [2,{2,3}], [3,{2,3}], [1,{1,2,3}], [2,{1,2,3}] and [3,{1,2,3}]
# beta for [1,{1}], [2,{2}] and [3,{3}] are omitted as obtained as [1,{1}] = 1 - [1,{1,2}]+[1,{1,3}]+[1,{1,2,3}] etc...

interior_alm_3d <- function(w,alpha,beta){

k <- c(1,2,3)

part1 <- (alpha[4]-1) * (2*alpha[4]-1)
part2 <- prod( beta[6+k]^alpha[4] * w[k]^(-alpha[4]-1) )
part3 <- sum( (beta[6+k]/w[k])^alpha[4] )^(1/alpha[4]-3)

return(part1*part2*part3)
}

edges_alm_3d <- function(w,alpha,beta,s,t){ # mass on the edge linking s-th and t-th components (in increasing order)

if(t<s){return('t cannot be less than s')}

if(s==1 && t==2){a=alpha[1];b1=beta[1];b2=beta[2]} ## s=1, t=2 or s=2, t=1
if(s==1 && t==3){a=alpha[2];b1=beta[3];b2=beta[4]} ## s=1, t=3 or s=3, t=1
if(s==2 && t==3){a=alpha[3];b1=beta[5];b2=beta[6]} ## s=3, t=2 or s=2, t=3
w1=w[s];w2=w[t];

part1 <- (a-1) * (b1*b2)^a * (w1*w2)^(-a-1)
part2 <- ( (b1/w1)^a + (b2/w2)^a )^(1/a-2)

return(part1*part2)
}

corners_alm_3d <- function(w,alpha,beta,s){ # mass on the corner of the s-th component
if(s==1){return(1-beta[1]-beta[2]-beta[4])}
if(s==2){return(1-beta[2]-beta[5]-beta[8])}
if(s==3){return(1-beta[4]-beta[6]-beta[9])}
}

dens_alm_3d <- function(w,alpha,beta,c){
if(length(alpha)!=4){return(stop("Wrong length of parameter alpha"))}
if(length(beta)!=9){return(stop("Wrong length of parameter beta"))}
if(beta[1]+beta[3]+beta[7]>1){return(1e-50)} # see conditions on the beta parameters
if(beta[2]+beta[5]+beta[8]>1){return(1e-50)}
if(beta[4]+beta[6]+beta[9]>1){return(1e-50)}
if(any(alpha < 1) || any(beta < 0) || any(beta > 1)){return(1e-50)}

if(c==0){return(interior_alm_3d(w,alpha,beta))}

if(sum(w<c) == 2){ # then we are in a corner
ind <- which(w > c)
return(corners_alm_3d(w,alpha,beta,ind)/c^2)
} else if(sum(w<c) == 1){
ind <- which(w >= c)
w2 <- w[ind]/sum(w[ind])
edge_surface <- c*sqrt(3)*(1-2*c)/2

if(w[1]<=1-c && w[2]<=1-c && w[3]<=c && w[1]>=(1-w[2])/2 && w[1]>=1-2*w[2]){ #EDGE {1,2}
edg01 <- integrate(Vectorize(function(x){edges_alm_3d(w=c(x,1-x,0), alpha=alpha, beta=beta, s=1, t=2)}), lower=0, upper=1)\$value
return(edges_alm_3d(c(w2,w[3]),alpha,beta,1,2) * edg01 / edge_surface)
}

if(w[1]<=1-c && w[2]<=c && w[3]<=1-c && w[2]<=w[1] && w[1]<=1-2*w[2]){ #EDGE {1,3}
edg01 <- integrate(Vectorize(function(x){edges_alm_3d(w=c(x,0,1-x), alpha=alpha, beta=beta, s=1, t=3)}), lower=0, upper=1)\$value
return(edges_alm_3d(c(w2[1],w[2],w2[2]),alpha,beta,1,3) * edg01 / edge_surface)
}

if(w[1]<=c && w[2]<=1-c && w[3]<=1-c && w[1]<=w[2] && w[1]<=(1-w[2])/2){ #EDGE {2,3}
edg01 <- integrate(Vectorize(function(x){edges_alm_3d(w=c(0,x,1-x), alpha=alpha, beta=beta, s=2, t=3)}), lower=0, upper=1)\$value
return(edges_alm_3d(c(w[1],w2),alpha,beta,2,3) * edg01 / edge_surface)
}
} else {
int01 <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_alm_3d(c(x,y,1-x-y),alpha=alpha, beta=beta)), lower=0, upper=1-y )\$value), lower=0, upper=1)\$value
intc <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_alm_3d(c(x,y,1-x-y),alpha=alpha, beta=beta)), lower=c, upper=1-y-c )\$value), lower=c, upper=1-2*c)\$value

return(interior_alm_3d(w,alpha,beta)*int01/intc)
}
}

### Angular density for the Asymmetric Logistic model on the 2 and 3 dimensional simplex

xvect = as.double(as.vector(t(x)))
if (is.vector(x)) {
dim = as.integer(length(x))
n = as.integer(1)
if(round(sum(x),7) !=1){ stop("Data is not angular") }
if(dim==2){ result <- dens_alm_2d(x,alpha,beta,c)}
if(dim==3){ result <- dens_alm_3d(x,alpha,beta,c)}
}
else {
dim = as.integer(ncol(x))
n = as.integer(nrow(x))
if (sum(apply(x,1,sum)) != n){ stop("Data is not angular") }
if (vectorial) {
result = double(n)
if(dim==2){ result <- apply(x,1,function(y){dens_alm_2d(y,alpha,beta,c)}) }
if(dim==3){ result <- apply(x,1,function(y){dens_alm_3d(y,alpha,beta,c)}) }
} else { # vectorial=FALSE mean we return the likelihood function
result = as.double(1)
if(dim==2){ result <- prod(apply(x,1,function(y){dens_alm_2d(y,alpha,beta,c)})) }
if(dim==3){ result <- prod(apply(x,1,function(y){dens_alm_3d(y,alpha,beta,c)})) }
}
}
if(log)
return(log(result))
else return(result)
}

### Functions needed for the Extremal-t model

dens_et <- function (x, rho, mu, c, log, vectorial){

# in the following functions:
#
# rho is a vector of size choose(dim,2) which mean choose 2 from dim
# mu is of size 1 (degree of freedom)

interior_et_d <- function(w,rho,mu){ # mass on the interior of the d-dim simplex
p<-length(w);
k=p-1;
w.tilde<-rep(0,k);
Sigma<-diag(k);

if(any(w==0) || any(w==1)){return(1e-300)}

for(i in 1:k){
w.tilde[i] <- ((w[i+1]/w[1])^(1/mu)-rho[i])*sqrt((mu+1)/(1-rho[i]^2));
for(j in i:k){
if(i==j){Sigma[i,j]=Sigma[j,i]=1}else{
Sigma[i,j]=Sigma[j,i]=(rho[sum(k:(k-i+1))+j-i]-rho[i]*rho[j])/sqrt((1-rho[i]^2)*(1-rho[j]^2))
}
}
}

if(sum(eigen(Sigma)\$values<0)>=1){return(1e-50)} #Check if matrix is positive definite

deriv = (w[-1]/w[1])^(1/mu-1)/mu*sqrt((mu+1)/(1-rho[1:k]^2))
if(k==1){
return(dest(x=w.tilde, scale=Sigma, df=mu+1)*w[1]^(-p-1)*prod(deriv) )
}else{
return(dmest(x=w.tilde, scale=Sigma, df=mu+1)*w[1]^(-p-1)*prod(deriv) )
}
}

# Mass on the other subsets of the 2d case:

corners_et_2d <- function(w,rho,mu,s){ # mass on the corner of the s-th component

part1 <- pt( -rho * sqrt((mu+1)/(1-rho^2)) ,df=mu+1)

return(part1)
}

dens_et_2d <- function(w,rho,mu,c){
if(length(rho)!=1){return(stop("Wrong length of parameter rho"))}
if( (abs(rho) > 1) || (mu<1) ){return(1e-50)}

if(sum(w<1-c) == 1){ # then we are in a corner
ind <- which(w > c)
return(corners_et_2d(w,rho,mu,ind))
} else {
return(interior_et_d(w,rho,mu))
}
}

# Mass on the other subsets of the 3d case:

corners_et_3d <- function(w,rho,mu,s){ # mass on the corner of the s-th component

if(s==1){rho1=rho[1];rho2=rho[2];rho3=rho[3]}
if(s==2){rho1=rho[1];rho2=rho[3];rho3=rho[2]}
if(s==3){rho1=rho[2];rho2=rho[3];rho3=rho[1]}

R <- (rho3 - rho1*rho2)/sqrt((1-rho1^2) * (1-rho2^2))
S <- matrix( c(1,R,R,1), ncol=2)
if(sum(eigen(S)\$values<0)>=1){return(1e-50)} # Check if matrix R1 is positive definite

a = -rho1 * sqrt((mu+1)/(1-rho1^2))
b = -rho2 * sqrt((mu+1)/(1-rho2^2))

return( pmest(x=c(a,b), scale=S, df=mu+1) )
}

# mass on the edge linking the s-th and t-th components (in inceasing order)

edges_et_3d <- function(w,rho,mu,s,t){

if(s==1 && t==2){rho1=rho[1];rho2=rho[2];rho3=rho[3]}
if(s==1 && t==3){rho1=rho[2];rho2=rho[1];rho3=rho[3]}
if(s==2 && t==3){rho1=rho[3];rho2=rho[1];rho3=rho[2]}
x=w[s];y=w[t];

A1 <- sqrt((mu+1)/(1-rho1^2)) * ( (y/x)^(1/mu) -rho1 )
A2 <- sqrt((mu+1)/(1-rho1^2)) * ( (x/y)^(1/mu) -rho1 )
B1 <- - rho2 * sqrt((mu+1)/(1-rho2^2))
B2 <- - rho3 * sqrt((mu+1)/(1-rho3^2))
d1 <- sqrt((mu+1)/(1-rho1^2)) / mu * y^(1/mu-1) /x^(1/mu)
d2 <- sqrt((mu+1)/(1-rho1^2)) / mu * x^(1/mu-1) /y^(1/mu)
R1 <- (rho3 - rho1*rho2)/sqrt((1-rho1^2) * (1-rho2^2))
R2 <- (rho2 - rho1*rho3)/sqrt((1-rho1^2) * (1-rho3^2))

S1 <- matrix(c(1,R1,R1,1),ncol=2)
S2 <- matrix(c(1,R2,R2,1),ncol=2)

if(sum(eigen(S1)\$values<0)>=1){return(1e-50)} # Check if matrix R1 is positive definite
if(sum(eigen(S2)\$values<0)>=1){return(1e-50)} # Check if matrix R2 is positive definite
if(any( abs(c(R1,R2)) >1 ) ){return(1e-50)}

part11 <- dt(x=A1,df=mu+1) * pt(sqrt((mu+2)/(mu+1+A1^2))*(B1-R1*A1)/sqrt(1-R1^2),df=mu+2) * d1 / x^2
part12 <- -1 - 1/mu + y * d1 * (mu+2) * A1 / (mu+1+A1^2)
part13 <- dt(x=A1,df=mu+1) * dt(x=sqrt((mu+2)/(mu+1+A1^2))*(B1-R1*A1)/sqrt(1-R1^2),df=mu+2) * d1^2 * y / x^2
part14 <- sqrt((mu+2)/(1-R1^2)) * (R1*(mu+1)+B1*A1) / (mu+1+A1^2)^(3/2)

part21 <- dt(x=A2,df=mu+1) * pt(sqrt((mu+2)/(mu+1+A2^2))*(B2-R2*A2)/sqrt(1-R2^2),df=mu+2) * d2 / y^2
part22 <- -1 - 1/mu + x * d2 * (mu+2) * A2 / (mu+1+A2^2)
part23 <- dt(x=A2,df=mu+1) * dt(x=sqrt((mu+2)/(mu+1+A2^2))*(B2-R2*A2)/sqrt(1-R2^2),df=mu+2) * d2^2 * x / y^2
part24 <- sqrt((mu+2)/(1-R2^2)) * (R2*(mu+1)+B2*A2) / (mu+1+A2^2)^(3/2)

return( -(x+y)^3 * (part11*part12 + part13*part14 + part21*part22 + part23*part24))

}

dens_et_3d <- function(w,rho,mu,c){

if(length(rho)!=3){return(stop("Wrong length of parameter rho"))}
if(any(abs(rho)>=1) || mu<=0 ){return(1e-50)}

if(c==0){return(interior_et_d(w,rho,mu))}

if(sum(w<c) == 2){ # then we are in a corner
ind <- which(w > c)
return(corners_et_3d(w,rho,mu,ind)/c^2)
}else if(sum(w<c) == 1){
ind <- which(w >= c)
w2 <- w[ind]/sum(w[ind])
edge_surface <- c*sqrt(3)*(1-2*c)/2

if(w[1]<=1-c && w[2]<=1-c && w[3]<=c && w[1]>=(1-w[2])/2 && w[1]>=1-2*w[2]){ #EDGE {1,2}
edg01 <- integrate(Vectorize(function(x){edges_et_3d(w=c(x,1-x,0), rho=rho, mu=mu, s=1, t=2)}), lower=0, upper=1)\$value
return(edges_et_3d(c(w2,w[3]),rho,mu,1,2) * edg01 / edge_surface)
}

if(w[1]<=1-c && w[2]<=c && w[3]<=1-c && w[2]<=w[1] && w[2]<=(1-w[1])/2){ #EDGE {1,3}
edg01 <- integrate(Vectorize(function(x){edges_et_3d(w=c(x,0,1-x), rho=rho, mu=mu, s=1, t=3)}), lower=0, upper=1)\$value
return(edges_et_3d(c(w2[1],w[2],w2[2]),rho,mu,1,3) * edg01 / edge_surface)
}

if(w[1]<=c && w[2]<=1-c && w[3]<=1-c && w[1]<=w[2] && w[1]<=(1-w[2])/2){ #EDGE {2,3}
edg01 <- integrate(Vectorize(function(x){edges_et_3d(w=c(0,x,1-x), rho=rho, mu=mu, s=2, t=3)}), lower=0, upper=1)\$value
return(edges_et_3d(c(w[1],w2),rho,mu,2,3) * edg01 / edge_surface)
}
}	else {
int01 <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_et_d(c(x,y,1-x-y),rho=rho, mu=mu)), lower=0, upper=1-y )\$value), lower=0, upper=1)\$value
intc <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_et_d(c(x,y,1-x-y),rho=rho, mu=mu)), lower=c, upper=1-y-c )\$value), lower=c, upper=1-2*c)\$value
return(interior_et_d(w,rho,mu)*int01/intc)
}
}

### Angular density for the Extremal-t model on the 2 and 3 dimensional simplex

xvect = as.double(as.vector(t(x)))
if (is.vector(x)) {
dim = as.integer(length(x))
n = as.integer(1)
if(round(sum(x),7) != 1){ stop("Data is not angular")}
if(dim==2){ result <- dens_et_2d(x,rho,mu,c)}
if(dim==3){ result <- dens_et_3d(x,rho,mu,c)}
}
else {
dim = as.integer(ncol(x))
n = as.integer(nrow(x))

if (sum(apply(x,1,sum)) != n){ stop("Data is not angular") }

if (vectorial) {
result = double(n)
if(dim==2){ result <- apply(x,1,function(y){dens_et_2d(y,rho,mu,c)}) }
if(dim==3){ result <- apply(x,1,function(y){dens_et_3d(y,rho,mu,c)}) }
} else { # vectorial=FALSE mean we return the likelihood function
result = as.double(1)
if(dim==2){ result <- prod(apply(x,1,function(y){dens_et_2d(y,rho,mu,c)})) }
if(dim==3){ result <- prod(apply(x,1,function(y){dens_et_3d(y,rho,mu,c)})) }
}
}
if(log){
if(any(result==0)){
result[which(result==0)] <- log(-1e+50)
result[which(result!=0)] <- log(result[which(result!=0)])
return(result)
}else{
return(log(result))
}
}else{
return(result)
}
}

### Functions needed for the Extremal Skew-t model

dens_est <- function (x, rho, alpha, mu, c, log, vectorial){

## Preliminary functions

Sigma <- function(rho){
p <- round(uniroot(function(x){length(rho)-choose(x,2)},lower=1, upper=10)\$root)
Sig <- diag(p)
Sig[lower.tri(Sig)] = rho
Sig[row(Sig) < col(Sig)] = t(Sig)[row(Sig) < col(Sig)]
return(Sig)
}

Sigma_j <-function(rho,j){
Sig <- Sigma(rho)
return(Sig[-j,-j] - Sig[-j,j] %*% t(Sig[j,-j]) )
}

s_j <- function(rho,j){
p <- round(uniroot(function(x){length(rho)-choose(x,2)},lower=1, upper=10)\$root)
k=p-1;

sigma_j <- Sigma_j(rho,j)
M <- diag(k)
diag(M) = sqrt(diag(sigma_j))
return( M )
}

Sigma_bar_j <- function(rho,j){
sigma_j <- Sigma_j(rho,j)
sj <- s_j(rho,j)
return( solve(sj) %*% t(sigma_j) %*% solve(sj) )
}

alpha_tilde <- function(alpha,j){
return(t(alpha[-j]))
}

alpha_bar_j <- function(rho,alpha,j){
Sig <- Sigma(rho)
sigma_j <- Sigma_j(rho,j)
Alpha_tilde <- alpha_tilde(alpha,j)

num <- alpha[j] + Sig[j,-j] %*% t(Alpha_tilde)
denom <- sqrt( 1 + Alpha_tilde %*% t(sigma_j) %*% t(Alpha_tilde)  )
return(num/denom)
}

alpha_star_j <- function(alpha,j){
Alpha_tilde <- alpha_tilde(alpha,j)
sj <- s_j(rho,j)
return( Alpha_tilde %*% sj )
}

tau_star_j <- function(rho,alpha,j){
Sig <- Sigma(rho)
Alpha_tilde <- alpha_tilde(alpha,j)
return( sqrt(mu+1) * (alpha[j] + Sig[-j,j] %*% t(Alpha_tilde) )	 )
}

nu_p <- function(rho,alpha,mu,j){
Alpha_bar <- alpha_bar_j(rho,alpha,j)
return( pt(sqrt(mu+1)*Alpha_bar,df=mu+1) * 2^(mu/2-1) * gamma((mu+1)/2) / sqrt(pi) )
}

x_bar <- function(x,rho,alpha,mu,j){
return( x * nu_p(rho,alpha,mu,j) )
}

comp <- function(z1,z2,z3,rho,alpha,mu,j,k){

#function corresponding to the j-th component of I_k
# Gives x1,y1, x2,y2, x3,y3
# a1,b1, a2,b2, a3,b3

#j = 1 or 2
#k = 1,2 or 3

if(j==1 & k==1){corr=rho[1];x=x_bar(z2,rho,alpha,mu,2); y=x_bar(z1,rho,alpha,mu,1)}
if(j==1 & k==2){corr=rho[1];x=x_bar(z1,rho,alpha,mu,1); y=x_bar(z2,rho,alpha,mu,2)}
if(j==1 & k==3){corr=rho[2];x=x_bar(z1,rho,alpha,mu,1); y=x_bar(z3,rho,alpha,mu,3)}

if(j==2 & k==1){corr=rho[2];x=x_bar(z3,rho,alpha,mu,3); y=x_bar(z1,rho,alpha,mu,1)}
if(j==2 & k==2){corr=rho[3];x=x_bar(z3,rho,alpha,mu,3); y=x_bar(z2,rho,alpha,mu,2)}
if(j==2 & k==3){corr=rho[3];x=x_bar(z2,rho,alpha,mu,2); y=x_bar(z3,rho,alpha,mu,3)}

if(z1==0 && z2==0 && z3==0){
return( -corr * sqrt((mu+1)/(1-corr^2)) )
}else{
return( sqrt((mu+1)/(1-corr^2))*((x/y)^(1/mu)-corr) )
}
}

R_j <- function(rho,alpha,j,matrix=TRUE){
sigma_bar <- Sigma_bar_j(rho,j)
d <- ncol(sigma_bar)
delta <- delta_j(rho,alpha,j)

S <- diag(d+1)
S[(1:d),(1:d)] <- sigma_bar
S[(1:d),d+1] <- S[d+1,(1:d)] <- -delta
if (matrix==TRUE){return(S)}else{
seq <- vector("numeric")
for(i in 1:d){
seq <- c(seq,S[(i+1):(d+1),i])
}
return(seq)
}
}

delta_j <- function(rho,alpha,j){

sigma_bar <- Sigma_bar_j(rho,j)
alpha_star <- alpha_star_j(alpha,j)
return( (alpha_star %*% sigma_bar )/ as.numeric(sqrt( 1 + alpha_star %*% sigma_bar %*% t(alpha_star) )) )
}

tau_bar_j <- function(rho,alpha,j){
tau_star <- tau_star_j(rho,alpha,j)
alpha_star <- alpha_star_j(alpha,j)
sigma_bar <- Sigma_bar_j(rho,j)
return( tau_star / sqrt( 1 + alpha_star %*% sigma_bar %*% t(alpha_star) ) )
}

# in the following functions:
#
# rho is a vector of size choose(dim,2) which mean choose 2 from dim (correlation coefficients)
# alpha is a vector of size dim (skewness parameters)
# mu is of size 1 (degree of freedom)

## Mass on the interior of the simplex

interior_skewt_d <- function(w,rho,alpha,mu){ # mass on the interior of the d-dim simplex
p<-length(w);
k=p-1;
w.tilde<-rep(0,k);

cond <- sapply(1:p,function(j){alpha_tilde(alpha,j) %*% t(Sigma_j(rho,j)) %*% t(alpha_tilde(alpha,j))})

if( any(cond < -1)){return(1e-50)}else{
sigma_bar1 <- Sigma_bar_j(rho,1)
alpha_star1 <- alpha_star_j(alpha,1)
tau_star1 <- tau_star_j(rho,alpha,1)

w_bar <- sapply(1:p,function(j){w[j]*nu_p(rho,alpha,mu,j)})

for(i in 1:k){
w.tilde[i] <- ((w_bar[i+1]/w_bar[1])^(1/mu)-rho[i])*sqrt((mu+1)/(1-rho[i]^2));
}

# if( alpha_star1 %*% sigma_bar1 %*% t(alpha_star1) < -1){return(1e-50)}
# if( t(w.tilde) %*% solve(sigma_bar1) %*% w.tilde < -1){return(1e-50)}

deriv = (w_bar[-1]/w_bar[1])^(1/mu-1) / mu * sqrt((mu+1)/(1-rho[1:k]^2)) * sapply(2:p,function(x){nu_p(rho,alpha,mu,x)}) / nu_p(rho,alpha,mu,1)

if(k==1){
return(dest(x=w.tilde, scale=sigma_bar1, shape=t(alpha_star1), extended=tau_star1, df=mu+1)*w[1]^(-p-1)*prod(deriv) )
}else{
return(dmest(x=w.tilde, location=rep(0,2), scale=sigma_bar1, shape=t(alpha_star1), extended=tau_star1, df=mu+1) * w[1]^(-p-1) * prod(deriv) )
}

}
}

## mass on the corner of the s-th component of the 2-d simplex

corners_skewt_2d <- function(w,rho,alpha,mu,s){

alpha_star <- alpha_star_j(alpha,s)
tau_star <- tau_star_j(rho,alpha,s)
part1 <- pest(-rho * sqrt((mu+1)/(1-rho^2)),shape=alpha_star, extended=tau_star ,df=mu+1)

return(part1)
}

## density on 2-d simplex

dens_skewt_2d <- function(w,rho,alpha,mu,c){
if(length(rho)!=1){return(stop("Wrong length of parameter rho"))}
if( (abs(rho) > 1) || (mu<1) ){return(1e-50)}

if(sum(w<c) == 1){ # then we are in a corner
ind <- which(w > c)
return(corners_skewt_2d(w,rho,alpha,mu,ind))
} else {
return(interior_skewt_d(w,rho,alpha,mu))
}
}

## mass on the corner of the s-th component of the 3-d simplex

corners_skewt_3d <- function(w,rho,alpha,mu,s){

s_bar <- Sigma_bar_j(rho,s)
al=t(alpha_star_j(alpha,s))
if(t(al) %*% s_bar %*% al < -1 ){return(1e-50)}

up <- c(comp(0,0,0,rho,alpha,mu,1,s), comp(0,0,0,rho,alpha,mu,2,s))
return(pmest(x=up, location=rep(0,2), scale=s_bar, shape=al, extended=tau_star_j(rho,alpha,s), df=mu+1 ))
}

## mass on the edge between the s-th and t-th components of the 3-d simplex

edges_skewt_3d <- function(w,rho,alpha,mu,s,t){

sigma_j1 <- Sigma_j(rho,s)
Alpha_tilde1 <- alpha_tilde(alpha,s)
sigma_j2 <- Sigma_j(rho,t)
Alpha_tilde2 <- alpha_tilde(alpha,t)
alpha_star1 <- alpha_star_j(alpha,s)
sigma_bar1 <- Sigma_bar_j(rho,s)
alpha_star2 <- alpha_star_j(alpha,t)
sigma_bar2 <- Sigma_bar_j(rho,t)
if( Alpha_tilde1 %*% t(sigma_j1) %*% t(Alpha_tilde1)< -1){return(1e-50)}
if( Alpha_tilde2 %*% t(sigma_j2) %*% t(Alpha_tilde2)< -1){return(1e-50)}
if(	alpha_star1 %*% sigma_bar1 %*% t(alpha_star1)< -1){return(1e-50)}
if(	alpha_star2 %*% sigma_bar2 %*% t(alpha_star2)< -1){return(1e-50)}

ind= c(1,2,3)
w[which((ind != s)  & (ind != t))]=0

if(s==1 && t==2){
a1_p <- comp(w[1],w[2],w[3],rho,alpha,mu,1,s)
b1_p <- comp(w[1],w[2],w[3],rho,alpha,mu,2,s)
R1 <- R_j(rho,alpha,s,FALSE)

a2_p <- comp(w[1],w[2],w[3],rho,alpha,mu,1,t)
b2_p <- comp(w[1],w[2],w[3],rho,alpha,mu,2,t)
R2 <- R_j(rho,alpha,t,FALSE)
}
if(s==1 && t==3){
a1_p <- comp(w[1],w[2],w[3],rho,alpha,mu,2,s)
b1_p <- comp(w[1],w[2],w[3],rho,alpha,mu,1,s)
R1 <- R_j(rho,alpha,s,FALSE)
r1 <- R1[2];
r2 <- R1[3];
R1[2] <- r2;
R1[3] <- r1;

a2_p <- comp(w[1],w[2],w[3],rho,alpha,mu,1,t)
b2_p <- comp(w[1],w[2],w[3],rho,alpha,mu,2,t)
R2 <- R_j(rho,alpha,t,FALSE)
}
if(s==2 && t==3){
a1_p <- comp(w[1],w[2],w[3],rho,alpha,mu,2,s)
b1_p <- comp(w[1],w[2],w[3],rho,alpha,mu,1,s)
R1 <- R_j(rho,alpha,s,FALSE)
r1 <- R1[2];
r2 <- R1[3];
R1[2] <- r2;
R1[3] <- r1;

a2_p <- comp(w[1],w[2],w[3],rho,alpha,mu,2,t)
b2_p <- comp(w[1],w[2],w[3],rho,alpha,mu,1,t)
R2 <- R_j(rho,alpha,t,FALSE)
rr1 <- R2[2];
rr2 <- R2[3];
R2[2] <- rr2;
R2[3] <- rr1;
}

if( R1[1]^2 >1 || R1[2]^2>1 || R2[1]^2 >1 || R2[2]^2>1 ){return(1e-50)}

c1 <- tau_bar_j(rho,alpha,s)
u1_p <- c(a1_p,b1_p,c1)

R1_p <- diag(2)
R1_p[1,2] <- R1_p[2,1] <- (R1[3]-R1[1]*R1[2])/sqrt((1-R1[1]^2)*(1-R1[2]^2))

c2 <- tau_bar_j(rho,alpha,t)
u2_p <- c(a2_p,b2_p,c2)

R2_p <- diag(2)
R2_p[1,2] <- R2_p[2,1] <- (R2[3]-R2[1]*R2[2])/sqrt((1-R2[1]^2)*(1-R2[2]^2))

x1_bar <- x_bar(w[s],rho,alpha,mu,s)
x2_bar <- x_bar(w[t],rho,alpha,mu,t)

v1_1 <- (b1_p-R1[1]*a1_p)/sqrt(1-R1[1]^2) * sqrt((mu+2)/(mu+1+a1_p^2))
v2_1 <- (c1-R1[2]*a1_p)/sqrt(1-R1[2]^2) * sqrt((mu+2)/(mu+1+a1_p^2))

v1_1_p <- (a1_p-R1[1]*b1_p)/sqrt(1-R1[1]^2) * sqrt((mu+2)/(mu+1+b1_p^2))

v1_2 <- (b2_p-R2[1]*a2_p)/sqrt(1-R2[1]^2) * sqrt((mu+2)/(mu+1+a2_p^2))
v2_2 <- (c2-R2[2]*a2_p)/sqrt(1-R2[2]^2) * sqrt((mu+2)/(mu+1+a2_p^2))

v1_2_p <- (a2_p-R2[1]*b2_p)/sqrt(1-R2[1]^2) * sqrt((mu+2)/(mu+1+b2_p^2))

if(s==1 & t==2){rho12=rho[1];rho13=rho[2];rho23=rho[3]}
if(s==1 & t==3){rho12=rho[2];rho13=rho[1];rho23=rho[3]}
if(s==2 & t==3){rho12=rho[3];rho13=rho[1];rho23=rho[2]}

part1 <- pt(c1,df=mu+1)
part11 <- -dt(a1_p,df=mu+1) * pmest(x=c( v1_1, v2_1), scale=R1_p, df=mu+2) * sqrt((mu+1)/(1-rho12^2)) * ( x2_bar/x1_bar )^(1/mu-1)
part12 <- nu_p(rho,alpha,mu,s)^2 * nu_p(rho,alpha,mu,t) / (mu*x1_bar^3) * ( 1+ 1/mu )

p1 <- dt(a1_p,df=mu+1)/(mu+1+a1_p^2)
p2 <- (mu+2) * a1_p * pmest(x=c( v1_1, v2_1), scale=R1_p, df=mu+2 )

p3 <- dt(v1_1,df=mu+2) * sqrt((mu+2)/(1-R1[1]^2)) * (b1_p*a1_p+R1[1]*(mu+1))/sqrt(mu+1+a1_p^2)
p4num <- sqrt(mu+3) * ( (c1 -R1[2]*a1_p)*(1-R1[1]^2) - (R1[3]-R1[1]*R1[2])*(b1_p-R1[1]*a1_p) )
p4denom <- ((1-R1[1]^2)*(mu+1+a1_p^2)+(b1_p-R1[1]*a1_p)^2) * ((1-R1[1]^2)*(1-R1[2]^2)-(R1[3]-R1[1]*R1[2])^2)
p4 <- pt(p4num/sqrt(p4denom),df=mu+2)

p5 <- dt(v2_1,df=mu+2) * sqrt((mu+2)/(1-R1[2]^2)) * (c1*a1_p+R1[2]*(mu+1))/sqrt(mu+1+a1_p^2)
p6num <- sqrt(mu+3) * ( (b1_p -R1[1]*a1_p)*(1-R1[2]^2) - (R1[3]-R1[1]*R1[2])*(c1-R1[2]*a1_p) )
p6denom <- ((1-R1[2]^2)*(mu+1+a1_p^2)+(c1-R1[2]*a1_p)^2) * ((1-R1[1]^2)*(1-R1[2]^2)-(R1[3]-R1[1]*R1[2])^2)
p6 <- pt(p6num/sqrt(p6denom),df=mu+2)

part13 <- (p1 * (p2+ p3*p4 + p5*p6)) * (mu+1)/(1-rho12^2) * ( x2_bar/x1_bar )^(2/mu-1)
part14 <- nu_p(rho,alpha,mu,s)^2 * nu_p(rho,alpha,mu,t) / (mu^2 * x1_bar^3)

part2 <- pt(c2,df=mu+1)
part21 <- -dt(a2_p,df=mu+1) * pmest(x=c( v1_2, v2_2), scale=R2_p, df=mu+2 ) * sqrt((mu+1)/(1-rho12^2)) * ( x1_bar/x2_bar )^(1/mu-1)
part22 <- nu_p(rho,alpha,mu,s) * nu_p(rho,alpha,mu,t)^2 / (mu*x2_bar^3) * ( 1+ 1/mu )

P1 <- dt(a2_p,df=mu+1)/(mu+1+a2_p^2)
P2 <- (mu+2) * a2_p * pmest(x=c( v1_2, v2_2), scale=R2_p, df=mu+2 )

P3 <- dt(v1_2,df=mu+2) * sqrt((mu+2)/(1-R2[1]^2)) * (b2_p*a2_p+R2[1]*(mu+1))/sqrt(mu+1+a2_p^2)
P4num <- sqrt(mu+3) * ( (c2 -R2[2]*a2_p)*(1-R2[1]^2) - (R2[3]-R2[1]*R2[2])*(b2_p-R2[1]*a2_p) )
P4denom <- ((1-R2[1]^2)*(mu+1+a2_p^2)+(b2_p-R2[1]*a2_p)^2) * ((1-R2[1]^2)*(1-R2[2]^2)-(R2[3]-R2[1]*R2[2])^2)
P4 <- pt(P4num/sqrt(P4denom),df=mu+2)

P5 <- dt(v2_2,df=mu+2) * sqrt((mu+2)/(1-R2[2]^2)) * (c2*a2_p+R2[2]*(mu+1))/sqrt(mu+1+a2_p^2)
P6num <- sqrt(mu+3) * ( (b2_p -R2[1]*a2_p)*(1-R2[2]^2) - (R2[3]-R2[1]*R2[2])*(c2-R2[2]*a2_p) )
P6denom <- ((1-R2[2]^2)*(mu+1+a2_p^2)+(c2-R2[2]*a2_p)^2) * ((1-R2[1]^2)*(1-R2[2]^2)-(R2[3]-R2[1]*R2[2])^2)
P6 <- pt(P6num/sqrt(P6denom),df=mu+2)

part23 <- (P1 * (P2+ P3*P4 + P5*P6)) * (mu+1)/(1-rho12^2) * ( x1_bar/x2_bar )^(2/mu-1)
part24 <- nu_p(rho,alpha,mu,s) * nu_p(rho,alpha,mu,t)^2 / (mu^2 * x2_bar^3)

return( -(w[s]+w[t])^3 * ((part11*part12+part13*part14)/part1 + (part21*part22+part23*part24)/part2))

}

## density on 3-d simplex

dens_skewt_3d <- function(w,rho,alpha,mu,c){

if(length(rho)!=3){return(stop("Wrong length of parameter rho"))}
if(length(alpha)!=3){return(stop("Wrong length of parameter rho"))}
if(any(abs(rho)>=1) || mu<=0){return(1e-50)}

if(c==0){return(interior_skewt_d(w,rho,alpha,mu))}

if(sum(w<c) == 2){ # then we are in a corner
ind <- which(w > c)
return(corners_skewt_3d(w,rho,alpha,mu,ind)/c^2)
} else
if(sum(w<c) == 1){
ind <- which(w >= c)
w2 <- w[ind]/sum(w[ind])
edge_surface <- c*sqrt(3)*(1-2*c)/2

if(w[1]<=1-c && w[2]<=1-c && w[3]<=c && w[1]>=(1-w[2])/2 && w[1]>=1-2*w[2]){ #EDGE {1,2}
edg01 <- integrate(Vectorize(function(x){edges_skewt_3d(w=c(x,1-x,0), rho=rho, alpha=alpha, mu=mu, s=1, t=2)}), lower=0, upper=1)\$value
return(edges_skewt_3d(c(w2,w[3]),rho,alpha,mu,1,2) * edg01 / edge_surface)
}

if(w[1]<=1-c && w[2]<=c && w[3]<=1-c && w[2]<=w[1] && w[1]<=1-2*w[2]){ #EDGE {1,3}
edg01 <- integrate(Vectorize(function(x){edges_skewt_3d(w=c(x,0,1-x), rho=rho, alpha=alpha, mu=mu, s=1, t=3)}), lower=0, upper=1)\$value
return(edges_skewt_3d(c(w2[1],w[2],w2[2]),rho,alpha,mu,1,3) * edg01 / edge_surface)
}

if(w[1]<=c && w[2]<=1-c && w[3]<=1-c && w[1]<=w[2] && w[1]<=(1-w[2])/2){ #EDGE {2,3}
edg01 <- integrate(Vectorize(function(x){edges_skewt_3d(w=c(0,x,1-x), rho=rho, alpha=alpha, mu=mu, s=2, t=3)}), lower=0, upper=1)\$value
return(edges_skewt_3d(c(w[1],w2),rho,alpha,mu,2,3) * edg01 / edge_surface)
}
}	else {
int01 <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_skewt_d(c(x,y,1-x-y),rho=rho, alpha=alpha, mu=mu)), lower=0, upper=1-y )\$value), lower=0, upper=1)\$value
intc <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_skewt_d(c(x,y,1-x-y),rho=rho, alpha=alpha, mu=mu)), lower=c, upper=1-y-c )\$value), lower=c, upper=1-2*c)\$value

return(interior_skewt_d(w,rho,alpha,mu)*int01/intc)
}
}

### Angular density for the Extremal Skew-t model on the 2 and 3 dimensional simplex

xvect = as.double(as.vector(t(x)))
if (is.vector(x)) {
dim = as.integer(length(x))
n = as.integer(1)
if(round(sum(x),7) != 1){ stop("Data is not angular")}
if(dim==2){ result <- dens_skewt_2d(x,rho,alpha,mu,c)}
if(dim==3){ result <- dens_skewt_3d(x,rho,alpha,mu,c)}
}
else {
dim = as.integer(ncol(x))
n = as.integer(nrow(x))

if (sum(apply(x,1,sum)) != n){ stop("Data is not angular") }

if (vectorial) {
result = double(n)
if(dim==2){ result <- apply(x,1,function(y){dens_skewt_2d(y,rho,alpha,mu,c)}) }
if(dim==3){ result <- apply(x,1,function(y){dens_skewt_3d(y,rho,alpha,mu,c)}) }
} else { # vectorial=FALSE mean we return the likelihood function
result = as.double(1)
if(dim==2){ result <- prod(apply(x,1,function(y){dens_skewt_2d(y,rho,alpha,mu,c)})) }
if(dim==3){ result <- prod(apply(x,1,function(y){dens_skewt_3d(y,rho,alpha,mu,c)})) }
}
}
if(log){
if(result==0){return(log(1e-50))}else{
return(log(result))
}
}else return(result)
}

# Dimension of the density
if (is.vector(x)) { d <- as.integer(length(x))}else{d <- as.integer(ncol(x)) }

if(model=='Pairwise'){
if(length(par)!= (choose(d,2)+1) ){stop('Wrong length of parameters')}
return(dens_pb(x=x, b=par[1:choose(d,2)], alpha=par[choose(d,2)+1], log=log,  vectorial=vectorial))
}
if(model=='Husler'){
if(length(par)!= choose(d,2) ){stop('Wrong length of parameters')}
return(dens_hr(x=x, lambda=par, log=log,  vectorial=vectorial))
}
if(model=='Dirichlet'){
if(length(par)!= d ){stop('Wrong length of parameters')}
return(dens_di(x=x, para=par, log=log,  vectorial=vectorial))
}
if(model=='Extremalt'){
if(length(par)!= (choose(d,2)+1) ){stop('Wrong length of parameters')}
if(is.null(c)){stop('c needs to be specified')}
return(dens_et(x=x, rho=par[1:choose(d,2)], mu=par[choose(d,2)+1], c=c, log=log,  vectorial=vectorial))
}
if(model=='Skewt'){
if(length(par)!= (choose(d,2)+d+1) ){stop('Wrong length of parameters')}
if(is.null(c)){stop('c needs to be specified')}
return(dens_est(x=x, rho=par[1:choose(d,2)], alpha=par[choose(d,2)+1:d],mu=par[choose(d,2)+d+1], c=c, log=log,  vectorial=vectorial))
}
if(model=='Asymmetric'){
if(is.null(c)){stop('c needs to be specified')}
if(d==2){
if(length(par)!=3 ){
stop('Wrong length of parameters')
}else{
return(dens_al(x=x, alpha=par[1], beta=par[2:3], c=c, log=log,  vectorial=vectorial))
}
}
if(d==3){
if(length(par)!=13 ){
stop('Wrong length of parameters')
}else{
return(dens_al(x=x, alpha=par[1:4], beta=par[5:13], c=c, log=log,  vectorial=vectorial))
}
}
}

}

### Definition of the angular density function and its rescaled version, i.e.
#
#  h*(w) := A''(w)/(A'(1) - A'(0))
#
dh <- function(w, beta, mixture = FALSE){
k <- length(beta) - 1
j <- 1:(k-1)
const <- 2/(k * (2-beta[2]-beta[k]))
res <- diff(diff(beta)) * dbeta(w, j, k-j)
if(mixture) return(k/2 * const * sum(res))
else return(k/2 * sum(res))
}

### Definition of the angular probability measure

ph <- function(w, beta)
{
k <- length(beta) - 1
j <- 1:k
#  res <- diff(beta) * dbeta(w, j, k-j+1)
res <- 0.5 * (diff(beta) + 1/k) * dbeta(w, j, k-j+1)
return(sum(res))
}

### Simulate from a mixture of beta densities
#
#
#
rh.mixt <- function(N=1000, beta){
k <- length(beta) - 1
# Sample N random uniforms U
U = runif(N)

#Variable to store the samples from the mixture distribution
rand.samples = rep(NA,N)

#probabilities mixture
prob <- diff(diff(beta)) / (2-beta[2]-beta[k])

#Sampling from the mixture
rand.samples <- sapply(1:N, function(i, prob, k, U){
j <- min(which((U[i]<cumsum(prob)))); return(rbeta(1,j,k-j))},
prob, k, U)

return(rand.samples)
}

### Simulate from an angular density
#
#
#
rh <- function(N=1000, beta){
k <- length(beta) - 1
#Sample N random uniforms U
U = runif(N)

p0 <- (1+k*(beta[2]-1))/2
p1 <- (1-k*(1-beta[k]))/2
prob <- c(p0, 1-p0-p1, p1)

rand.samples <- sapply(1:N, function(i, prob, U){
j <- min(which((U[i]<cumsum(prob))));
return(switch(j, 0, rh.mixt(1,beta),1))},
prob,U)
}

### Estimation and generation from angular density

angular <- function(data, model, n, dep, asy, alpha, beta, df, seed, k, nsim, plot=TRUE, nw=100){

w <- seq(0.00001, .99999, length=nw)

# Simulation of synthetic data
if(missing(data)){
if(missing(model) || missing(n) || missing(dep) || missing(seed) ){stop("model, n. dep and missing must be specified when data is not provided")}
# warning("Data not provided and generated according to the parameters provided")

if(any(model==c("log", "alog", "hr", "neglog", "aneglog", "bilog", "negbilog", "ct", "amix"))){
set.seed(seed)
data <- rbvevd(n=n, dep = dep, asy, alpha, beta, model = model, mar1 = c(1,1,1))
Atrue <- abvevd(w, dep=dep, asy, alpha, beta, model=model) # True pickands dependence function
htrue <- hbvevd(1-w, dep=dep, asy, alpha, beta, model=model, half=TRUE) # True angular density
if(any(model==c("alog","aneglog"))){Htrue <- (1-asy)/2}
if(model=="amix"){Htrue <- c(1-alpha-beta, 1-alpha-2*beta)/2}
}else if(model=="Extremalt"){
set.seed(seed)
data <- r_extr_mod(model, n=n, param=c(dep, df))
Atrue <- rep(NA,nw)
for(i in 1:nw){Atrue[i] <- pk.extst(c(w[i],1-w[i]), c(dep,0,0,df))}
htrue <- dens(x=cbind(w,1-w), model=model, par=c(dep,df), c=0, log=FALSE, vectorial=TRUE)/2
Htrue <- dens(x=cbind(c(1,0),c(0,1)), model=model, par=c(dep,df), c=0, log=FALSE, vectorial=TRUE)/2
}
}else{
if(!is.matrix(data) || ncol(data)!=2){stop("data must be a matrix with 2 columns")}
n <- nrow(data)
model <- dep <- seed <- NULL
Atrue <- htrue <- 0
warning("True Pickands function and angular density not computed as data is provided and true model is unsure")
}

if(missing(k)){stop("k, the polynomial order must be specified")}
if(missing(nsim)){stop("nsim, the number of of generation from the estimated angular density must be specified")}

# Compute the angles:
wdata <- data[,1]/rowSums(data)

# Estimate the pickands function:
Aest <- beed(data, cbind(w, 1-w), 2, 'md', 'emp', k=k, plot=FALSE)
beta <- Aest\$beta

# Compute the angular density in the interior
hest <- sapply(1:nw, function(i, w, beta) dh(w[i], beta), w, beta)

# Compute the angular measure
Hest <- sapply(1:nw, function(i, w, beta) ph(w[i], beta), w, beta)

# Compute the point masses on the corners
p0 <- (1+k*(beta[2]-1))/2
p1 <- (1-k*(1-beta[k]))/2

# simulate from the semiparametric angular density
wsim <- rh(nsim, beta)

if(plot){
par(mai = c(0.85, 0.85, 0.1, 0.1), mgp = c(2.8, 1, 0), cex.axis=2, cex.lab=2)
hist(wsim[wsim!=0 & wsim!=1], freq=FALSE, ylim=c(0,3.5), xlim=c(0,1), xlab='w', ylab='h(w)', main="",col="gray")
lines(w, hest, col=1, lty=2, lwd=2.5)
if(is.vector(htrue)){ lines(w[2:99], htrue[2:99],col=2, lty=1, lwd=1.5)}
points(1,sum(wsim==1)/nsim, pch=19, cex=2)
points(0,sum(wsim==0)/nsim, pch=19, cex=2)
if(any(model==c('alog','aneglog','amix','Extremalt'))){
points(0, Htrue[1] , pch=21, cex=2,col=2, lwd=2)
points(1, Htrue[2], pch=21, cex=2,col=2, lwd=2)
}
}

return(list(model=model, n=n, dep=dep, data=data, Aest=Aest, hest=hest, Hest=Hest, p0=p0, p1=p1, wsim=wsim, Atrue=Atrue, htrue=htrue))

}
```

## Try the ExtremalDep package in your browser

Any scripts or data that you put into this service are public.

ExtremalDep documentation built on Aug. 29, 2019, 9:03 a.m.