Nothing
#############################################################
### Authors: Giulia Marcon and Simone Padoan ###
### Emails: giulia.marcon@phd.unibocconi.it, ###
### simone.padoan@unibocconi.it ###
### Institution: Department of Decision Sciences, ###
### University Bocconi of Milan ###
### File name: ConfidenceBands.r ###
### Description: ###
### This file contains routines that compute ###
### (1-a)% bootstrap confident bands for the estimated ###
### Pickands dependence function ###
### Last change: 15/08/2016 ###
#############################################################
# Subroutine that computes resampled estimates
beed.boot <- function(data, x, d=3, est = c("ht", "md", "cfg", "pick"),
margin=c("emp", "est", "exp", "frechet", "gumbel"),
k = 13, nboot = 500, y = NULL,
print = FALSE){
nx <- nrow(x)
ndata <- nrow(data)
ddata <- ncol(data)
if(d!=ddata)
stop("`data' must be a matrix/data.frame with `d' columns")
bootsamp <- matrix(0, nrow=ndata, ncol=d)
if(is.null(y))
fit <- beed(data = data, x = x, d = d, est = est, k = k,
plot = FALSE, margin = margin)
else fit <- y
A.tilde <- fit$A
beta.tilde <- fit$beta
beta <- bootA <- NULL
for(i in 1 : nboot)
{
indx <- sample(1:ndata, ndata, replace = TRUE)
bootsamp <- data[indx, ]
mod <- beed(data = bootsamp, x = x, d = d, est = est, k = k,
plot = FALSE, margin = margin)
b <- mod$beta
beta <- cbind(beta,b)
bootA <- cbind(bootA, mod$A)
if(print){
print.i <- seq(0, nboot, by=100)
if(i %in% print.i) message(paste(c('iteration', i, 'out of', nboot,'\n')))
}
}
out=list(A=A.tilde,bootA=bootA,beta=beta)
invisible(out)
}
# Main routine that computes the confidence bands
beed.confband <- function(data, x, d=3, est = c("ht", "md", "cfg", "pick"),
margin=c("emp", "est", "exp", "frechet", "gumbel"),
k = 13, nboot = 500, y = NULL, conf = 0.95,
plot = FALSE, print = FALSE){
nx <- nrow(x)
ndata <- nrow(data)
ddata <- ncol(data)
if(d!=ddata)
stop("`data' must be a matrix/data.frame with `d' columns")
# Boostrap
if(is.null(y))
fit <- beed.boot(data = data, x = x, d = d, est = est, k = k,
nboot = nboot, print = print,
margin = margin)
else
fit <- beed.boot(data = data, x = x, d = d, est = est, k = k,
y = y, nboot = nboot, print = print, margin = margin)
A.tilde <- fit$A
bootA <- fit$bootA
beta <- fit$beta
alpha <- 1-conf
# Confidence bands pointwise on x
A.low.pointwise <- A.up.pointwise <- numeric(nx)
for (i in 1:nx){
ord <- sort(bootA[i,])
A.low.pointwise[i] <- ord[round(nboot*(alpha/2))]
A.up.pointwise[i] <- ord[round(nboot*(1-alpha/2))]
}
# Nnumber of coefficients
p <- nrow(beta)
# Confidence bands on each beta
low.beta <- up.beta <- numeric(p)
for (i in 1:p){
ord <- sort(beta[i,])
low.beta[i] <- ord[round(nboot*(alpha/2))]
up.beta[i] <- ord[round(nboot*(1-alpha/2))]
}
A.up.beta <- beed(data = data, x = x, d = d, est = est, k = k,
beta = up.beta, margin = margin)$A
A.low.beta <- beed(data = data, x = x, d = d, est = est, k = k,
beta = low.beta, margin = margin)$A
if (plot == TRUE){
if(d == 2){
plot(x[,1],x[,2],type='n',xlab='w',ylab='A(w)',ylim=c(.5,1))
polygon(c(0, 0.5, 1), c(1, 0.5, 1), lty = 1, lwd = 1, border = 'grey')
lines(x[,1],A.tilde,lty=1,col=1)
lines(x[,1],A.low.beta,lty=2,col=1)
lines(x[,1],A.up.beta,lty=2,col=1)
}
if(d == 3){
numg <- sqrt(nx)
xy <- seq(0,1,length=numg)
mat <- matrix(A.tilde, numg, numg)
matup <- matrix(A.up.beta, numg, numg)
matlow <- matrix(A.low.beta, numg, numg)
plot(xy, xy, type='n', xlab=expression(w[1]), ylab=expression(w[2]))
image(x=xy, y=xy, z=mat, col=heat.colors(numg),add=TRUE)
contour(x=xy, y=xy, z=mat, add=T,col='black',labcex=.6,lty=1)
contour(x=xy, y=xy, z=matup, add=T,col='black',labcex=.6,lty=2)
contour(x=xy, y=xy, z=matlow, add=T,col='black',labcex=.6,lty=2)
}
if (d >= 4)
warning("cannot plot in high dimensions")
}
out=list(A=A.tilde,A.up.beta=A.up.beta,A.low.beta=A.low.beta,
bootA=bootA,A.up.pointwise=A.up.pointwise,
A.low.pointwise=A.low.pointwise,up.beta=up.beta,
low.beta=low.beta)
invisible(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.