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

#### Documented in beed

```#######################################################
### Authors: Giulia Marcon and Simone Padoan        ###
### Emails: [email protected],        ###
### [email protected]                     ###
### Institution: Department of Decision Sciences,   ###
### University Bocconi of Milan                     ###
### File name: beed.r                               ###
### Description:                                    ###
### This file provides an estimate of the uni- and  ###
### multi-variate Pickands dependence function,     ###
### based on an Bernstein polynomials approximation ###
### as proposed in Marcon et al. (2016)             ###
### Last change: 15/08/2016                         ###
#######################################################

############################################################################
# Bernstein Estimation of Extremal Dependence.                           ###
# Pickands Dependence Function Estimate is provided.                     ###
# NB. This relies on the representation of the                           ###
# exponent measure function assuming exponential marginals if using the  ###
# Madogram, whereas use frechet marginals when using the other           ###
# estimators for the pilot estimates                                     ###
############################################################################
################################################################################
# INPUT:                                                                     ###
# x is the d-dimensional simplex                                             ###
# y is the pilot estimate. If missing, it is estimated                       ###
# est = c('md','cfg','ht','pick') choses the method for the initial estimate ###
# k is the order of polynomials                                              ###
# d is the dimension of data                                                 ###
# beta is the vector of coefficient for the Bernstei representation.         ###
#   If it's missing, it is estimated                                         ###
################################################################################

beed <- function(data, x, d=3, est = c("ht","cfg","md","pick"),
margin=c("emp", "est", "exp", "frechet", "gumbel"),
k = 13, y = NULL, beta = NULL,
plot = FALSE){

xdim <- dim(x)

# Check conditions
stop("`data' must be a matrix with `d' columns")
if (!is.matrix(x) || xdim[2] != d)
stop("`x' must be a matrix with `d' columns")
if (any(x < 0, na.rm = TRUE))
stop("`x' must be non-negative")
rs <- rowSums(x)
nx <- xdim[1]

if (any(rs < 0, na.rm = TRUE))
warning("row(s) of `x' which do not lay in the simplex are removed")

# Subroutine that provides index beta coefficients
# p(k,d) = choose(k+d-1,d-1) number of coefficients
index <- function(k, d){
beta.index <- expand.grid(rep(list(0:k),d))
beta.index.sum <- rowSums(beta.index)
restr <- which(beta.index.sum<=k)
v <- as.matrix(beta.index[restr,ncol(beta.index):1])
return(v)
}

# Bivariate Case
if(d==2){
# Create design matrix:
xx <- as.matrix(x[,-d])
vb <- index(k=k,d-1)
q <- nrow(vb) # q = choose(k+d-1,d-1)
Z <- bp(x=x,k=k,v=vb)

# Preliminary Estimation if missing
if(is.null(y)){

if(length(est)>1){
est='md'
warning("invalid argument for `est', madogram by default")
}
if(length(margin)>1){
margin='emp'
warning("invalid argument for `margin', empirical transformation dy default")
}

if(margin=='emp')
y <- switch(est, ht = abvnonpar(x=xx, data=data, method="pickands", d=d, madj=2, epmar = TRUE),
md = madogram(w = x, data = data, margin = "emp"),
pick = An.biv(x=data, w=xx, estimator = "Pickands"),
cfg = An.biv(x=data, w=xx, estimator = "CFG")
)
if(margin=='est')
y <- switch(est, ht = abvnonpar(x=xx, data=data, method="pickands", d=d, madj=2, epmar = FALSE),
md = madogram(w = x, data = data, margin = "est"),
pick = An.biv(x=data, w=xx, estimator = "Pickands"),
cfg = An.biv(x=data, w=xx, estimator = "CFG")
)
if((margin!="emp")&(margin!="est")){
y <- switch(est, ht = abvnonpar(x=xx, data=data, method="pickands", d=d, madj=2, epmar = FALSE),
md = madogram(w = x, data = data, margin = margin),
pick = An.biv(x=data, w=xx, estimator = "Pickands"),
cfg = An.biv(x=data, w=xx, estimator = "CFG")
)
if((est=="ht")|(est=="pick")|(est=="cfg"))
warning("Marginal distributions estimated")
}
}

n <- length(y)

# Estimation procedure, else assign beta vector input
if(is.null(beta)){

# Create the A, b0 matrix for constraints
# R1) Convexity
R1 <- matrix(0,k-1,k+1)
for (i in 1:(k-1))
R1[i,i:(i+2)] <- c(1,-2,1)
r1 <- rep(0,nrow(R1))

# R2) upper bound
# A(ei) = 1
R2 <- matrix(0,2*d,q)
vertix <- which(apply(vb, 1, function(x) any(x==k)))
R2[1,1] <- 1
R2[(d+1),1]<- -1

for(i in 1:(d-1)) {
R2[i+1,vertix[i]] <- 1
R2[i+1+d,vertix[i]] <- -1
}
r2 <- c(rep(1,d),rep(-1,d))

# R3) lower bound
# A(w) >= max(w)
nr <- d*(d-1)
R3 <- matrix(0,nr,q)

vertix_10 <- which(vb==1)
vertix_k_10 <- which(vb==k-1)
vertix <- c(vertix_10,vertix_k_10)

for(i in 1:nr)
R3[i,vertix[i]] <- 1
r3 <- rep(1-1/k,nr)

R = rbind(R1, R2, R3)
r0 = c(r1, r2, r3)

Dmat  <- t(Z)%*%Z
dvec  <- t(y)%*%Z
Rmat  <- t(R)

fit <- solve.QP(Dmat,dvec,Rmat,r0)
beta.tilde <- as.matrix(fit\$solution)

}
else beta.tilde <- beta

y.tilde <- Z%*%beta.tilde
yy <- y

if (plot == TRUE){
plot(x[,1],x[,2],type='n',xlab='t',ylab='A(t)',ylim=c(.5,1))
polygon(c(0, 0.5, 1), c(1, 0.5, 1), lty = 1, lwd = 1, border = 'grey')
lines(x[,1],y.tilde,lty=1,col=1)
}
}

# mulivariate case
else{
zero1 <- which(rs>1)
zero0 <- which(rs==0)
zeroNa <- which(is.na(rs))
if((length(zero1)==0)&(length(zeroNa)==0)&(length(zero0)==0))
X <- x
else{
X <- x[-c(zero1,zeroNa,zero0),]
}

# Create design matrix:
xx <- as.matrix(X[,-d])
vb <- index(k,d-1)
q <- nrow(vb)
Z <- bp(x=X,k=k,v=vb)
# Preliminary Estimation if missing
if(is.null(y)){

if(length(est)>1){
est='md'
warning("invalid argument for `est', madogram by default")
}
if(length(margin)>1){
margin='emp'
warning("invalid argument for `margin', empirical transformation dy default")
}

if(margin=='emp')
y <- switch(est, ht = amvnonpar(x=X, data=data, d=d, madj=2, epmar = TRUE),
md = madogram(w = X, data = data, margin="emp"),
cfg = An(x = data, w = X)\$CFG,
pick = An(x = data, w = X)\$P
)
if(margin=='est')
y <- switch(est, ht = amvnonpar(x=X, data = data, d=d, madj=2, epmar = FALSE),
md = madogram(w = X, data = data, margin = "est"),
cfg = An(x = data, w = X)\$CFG,
pick = An(x = data, w = X)\$P
)
if((margin!="emp")&(margin!="est")){
y <- switch(est, ht = amvnonpar(x=X, data = data, d=d, madj=2, epmar = FALSE),
md = madogram(w = X, data = data, margin = margin),
cfg = An(x = data, w = X)\$CFG,
pick = An(x = data, w = X)\$P
)
if((est=="ht")|(est=="pick")|(est=="cfg"))
warning("Marginal distributions estimated")
}
}

n <- length(y)

yy <- rep(NA,nx)
if((length(zero1)==0)&(length(zeroNa)==0)&(length(zero0)==0)) yy <- y
else yy[-c(zero1,zeroNa,zero0)] <- y

# Estimation procedure, else assign beta vector input
if(is.null(beta)){

# Create the A, b0 matrix for constraints
# R1) Convexity
#R1 <- convexity(d = d, k = k)
#r1 <- rep(0,nrow(R1))
R1 <- convexity_old(v = vb, d = d)
r1 <- rep(0,nrow(R1))

# R2) upper bound
# A(ei) = 1
R2 <- matrix(0,2*d,q)
vertix <- which(apply(vb, 1, function(x) any(x==k)))
R2[1,1] <- 1
R2[(d+1),1]<- -1

for(i in 1:(d-1)) {
R2[i+1,vertix[i]] <- 1
R2[i+1+d,vertix[i]] <- -1
}
r2 <- c(rep(1,d),rep(-1,d))

# R3) lower bound
# A(w) >= max(w)
if(d==3){
nr <- d*(d-1)
R3 <- matrix(0,nr,q)

v_10 <- which(apply(vb, 1, function(x) any(x==1) &
any(x==0) & sum(x)==1))
v_k_10 <- which(apply(vb, 1, function(x) any(x==k-1) &
any(x==0)& sum(x)==k-1))
v_k_11 <- which(apply(vb, 1, function(x) any(x==k-1) &
any(x==1)& sum(x)==k))
vertix <- c(v_10,v_k_10,v_k_11)

for(i in 1:nr)
R3[i,vertix[i]] <- 1
r3 <- rep(1-1/k,nr)
}
if(d>3){
R3 <- Z
r3 <- apply(X,1,max)
}

R = rbind(R1, R2, R3)
r0 = c(r1, r2, r3)

Dmat  <- t(Z)%*%Z
dvec  <- t(y)%*%Z
Rmat  <- t(R)

fit <- solve.QP(Dmat,dvec,Rmat,r0)
beta.tilde <- as.matrix(fit\$solution)

}
else beta.tilde <- beta

y.tilde <- rep(NA,nx)
if((length(zero1)==0)&(length(zeroNa)==0)&(length(zero0)==0)) y.tilde <- Z%*%beta.tilde
else y.tilde[-c(zero1,zeroNa,zero0)] <- Z%*%beta.tilde

if (plot == TRUE){
if(d == 3){
if((length(zero1)==0)&(length(zeroNa)==0)&(length(zero0)==0)){
numg <- round((-1 + sqrt(1 + 8*xdim[1]))/2)
mat <- matrix(NA, numg, numg)
mat[lower.tri(mat, diag = T)] <- y.tilde
}
else{
numg <- sqrt(nx)
mat <- matrix(y.tilde, numg, numg)
}
xy <- seq(0,1,length=numg)
plot(xy, xy, type='n', xlab=expression(t[1]), ylab=expression(t[2]))
}
if (d >= 4)
warning("can not plot in high dimensions")
}

}

med <- matrix(rep(1/d,d),ncol=d)
Zmed <- bp(med,v=vb,k=k)
extind <- as.numeric( d * (Zmed%*%beta.tilde) )

out <- list(beta=beta.tilde,A=y.tilde, Anonconvex=yy, extind=extind)
return(out)
}

########################################################################
# Define a multivariate Bernstein Polynomials on the Simplex         ###
# X %*% beta = Bernstein Representation                              ###
# Input: x is the d-dimensional simplex, v index set, k polynomial   ###
# order                                                              ###
# bp returns a design matrix X                                       ###
########################################################################

bp <- function(x,k,v=NULL)
{

x <- as.data.frame(x)
dd <- ncol(x)-1
xx <- x[,(1:dd)]
if (dd == 1){
#x <- as.matrix(sort(xx))
q <- nrow(x) # number of points of the grid
basis <- array(0, dim=c(q, k+1))

for(j in 0:k){
basis[,j+1] <- choose(k,j) * x[,1]^j * x[,2]^(k-j)}

return(basis)
}

else x <- as.matrix(xx)
q <- nrow(v)
n <- nrow(x)
b <- matrix(NA,n,q)
xv <- matrix(NA,n,dd)

for (i in 1:q) {
for (j in 1:dd) {
xv[,j] <- x[,j]^v[i,j]
}
b[,i] <- factorial(k)/( prod(factorial(v[i,])) * factorial(k-sum(v[i,]))) *
apply(xv,1,prod) * (1-rowSums(x))^(k-sum(v[i,]))
}
return(b)
}

########################################################################
########################## BIVARIATE CASE ##############################
bp2d <- function(w, k)
{
#ind <- index(k) # compute the index
q <- nrow(w) # number of points of the grid
basis <- array(0, dim=c(q, k+1))

for(j in 0:k){
basis[,j+1] <- choose(k,j) * w[,1]^j * w[,2]^(k-j)}

return(basis)
}

########################## BIVARIATE CASE ##############################
########################################################################

########################################################################
# Define a multivariate Unit-Simplex                                 ###
# Input: n is the length of the one-dimension simplex                ###
# simplex returns a matrix of the d-dimensional Simplex              ###                           ###
########################################################################

simplex <- function(d, n=50, a=0, b=1){
w <- seq(a, b, length = n)
if(d==2) return(cbind(w,1-w))
else{
ww <- expand.grid(rep(list(w), d-1))
os <- which(rowSums(ww)>1)
x <- cbind(ww, 1 - rowSums(ww))
x <- x[-os,]
colnames(x) = paste("x", 1:d, sep="")
return(as.matrix(x))
}
}
```

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