Multivariate_Extremes.R
In evd: Functions for Extreme Value Distributions

### R code from vignette source 'Multivariate_Extremes.Rnw'

###################################################
### code chunk number 1: rawdata
###################################################
options(show.signif.stars=FALSE)
library(evd); nn <- nrow(lossalae)
loss <- lossalae/1e+05; lts <- c(1e-04, 100)
plot(loss, log = "xy", xlim = lts, ylim = lts)


###################################################
### code chunk number 2: Multivariate_Extremes.Rnw:42-44
###################################################
ula <- apply(loss, 2, rank)/(nn + 1)
plot(ula)


###################################################
### code chunk number 3: Multivariate_Extremes.Rnw:49-50
###################################################
options(show.signif.stars=FALSE)
library(evd); nn <- nrow(lossalae)
loss <- lossalae/1e+05; lts <- c(1e-04, 100)
plot(loss, log = "xy", xlim = lts, ylim = lts)


###################################################
### code chunk number 4: asylogdfn
###################################################
abvevd(dep = 0.5, asy = c(1,1), model = "alog", plot = TRUE)
abvevd(dep = 0.5, asy = c(0.6,0.9), model = "alog", add = TRUE, lty = 2)
abvevd(dep = 0.5, asy = c(0.8,0.5), model = "alog", add = TRUE, lty = 3)


###################################################
### code chunk number 5: Multivariate_Extremes.Rnw:88-91
###################################################
abvevd(dep = -1/(-2), model = "neglog", plot = TRUE)
abvevd(dep = -1/(-1), model = "neglog", add = TRUE, lty = 2)
abvevd(dep = -1/(-0.5), model = "neglog", add = TRUE, lty = 3)


###################################################
### code chunk number 6: Multivariate_Extremes.Rnw:94-97
###################################################
abvevd(alpha = 1, beta = -0.2, model = "amix", plot = TRUE)
abvevd(alpha = 0.6, beta = 0.1, model = "amix", add = TRUE, lty = 2)
abvevd(alpha = 0.2, beta = 0.2, model = "amix", add = TRUE, lty = 3)


###################################################
### code chunk number 7: Multivariate_Extremes.Rnw:100-103
###################################################
abvevd(dep = 1/1.25, model = "hr", plot = TRUE)
abvevd(dep = 1/0.83, model = "hr", add = TRUE, lty = 2)
abvevd(dep = 1/0.5, model = "hr", add = TRUE, lty = 3)


###################################################
### code chunk number 8: Multivariate_Extremes.Rnw:108-109
###################################################
abvevd(dep = 0.5, asy = c(1,1), model = "alog", plot = TRUE)
abvevd(dep = 0.5, asy = c(0.6,0.9), model = "alog", add = TRUE, lty = 2)
abvevd(dep = 0.5, asy = c(0.8,0.5), model = "alog", add = TRUE, lty = 3)


###################################################
### code chunk number 9: Multivariate_Extremes.Rnw:121-129
###################################################
set.seed(131); cml <- loss[sample(nn),]
xx <- rep(1:50, each = 30); lts <- c(1e-04, 100)
cml <- cbind(tapply(cml[,1], xx, max), tapply(cml[,2], xx, max))
colnames(cml) <- colnames(loss)
plot(loss, log = "xy", xlim = lts, ylim = lts, col = "grey")
points(cml)
ecml <- -log(apply(cml,2,rank)/51)
plot(ecml)


###################################################
### code chunk number 10: nonpardfn
###################################################
pp <- "pickands"; cc <- "cfg"
abvnonpar(data = cml, epmar = TRUE, method = pp, plot = TRUE, lty = 3)
abvnonpar(data = cml, epmar = TRUE, method = pp, add = TRUE, madj = 1, lty = 2)
abvnonpar(data = cml, epmar = TRUE, method = pp, add = TRUE, madj = 2, lty = 4)
abvnonpar(data = cml, epmar = TRUE, method = cc, add = TRUE, lty = 1)


###################################################
### code chunk number 11: Multivariate_Extremes.Rnw:142-148
###################################################
m1 <- fbvevd(cml, asy1 = 1, model = "alog")
m2 <- fbvevd(cml, model = "log")
m3 <- fbvevd(cml, model = "bilog")
plot(m1, which = 4, nplty = 3)
plot(m2, which = 4, nplty = 3, lty = 2, add = TRUE)
plot(m3, which = 4, nplty = 3, lty = 4, add = TRUE)


###################################################
### code chunk number 12: Multivariate_Extremes.Rnw:153-154
###################################################
pp <- "pickands"; cc <- "cfg"
abvnonpar(data = cml, epmar = TRUE, method = pp, plot = TRUE, lty = 3)
abvnonpar(data = cml, epmar = TRUE, method = pp, add = TRUE, madj = 1, lty = 2)
abvnonpar(data = cml, epmar = TRUE, method = pp, add = TRUE, madj = 2, lty = 4)
abvnonpar(data = cml, epmar = TRUE, method = cc, add = TRUE, lty = 1)


###################################################
### code chunk number 13: Multivariate_Extremes.Rnw:164-167
###################################################
round(rbind(fitted(m2), std.errors(m2)), 3)
anova(m3, m2)
evind.test(cml, method = "score")


###################################################
### code chunk number 14: nonparqc
###################################################
lts <- c(0.01,100)
plot(loss, log = "xy", col = "grey", xlim = lts, ylim = lts)
points(cml); pp <- c(0.98,0.99,0.995)
qcbvnonpar(pp, data = cml, epmar = TRUE, mint = 30, add = TRUE)


###################################################
### code chunk number 15: Multivariate_Extremes.Rnw:185-186
###################################################
lts <- c(0.01,100)
plot(loss, log = "xy", col = "grey", xlim = lts, ylim = lts)
points(cml); pp <- c(0.98,0.99,0.995)
qcbvnonpar(pp, data = cml, epmar = TRUE, mint = 30, add = TRUE)


###################################################
### code chunk number 16: Multivariate_Extremes.Rnw:202-204 (eval = FALSE)
###################################################
## k0 <- bvtcplot(loss)$k0
## bvtcplot(loss, spectral = TRUE)


###################################################
### code chunk number 17: bvtc
###################################################
k0 <- bvtcplot(loss)$k0


###################################################
### code chunk number 18: Multivariate_Extremes.Rnw:213-214
###################################################
k0 <- bvtcplot(loss)$k0


###################################################
### code chunk number 19: Multivariate_Extremes.Rnw:224-228
###################################################
thresh <- apply(loss, 2, sort, decreasing = TRUE)[(k0+5)/2,]
mar1 <- fitted(fpot(loss[,1], thresh[1]))
mar2 <- fitted(fpot(loss[,2], thresh[2]))
rbind(mar1,mar2)


###################################################
### code chunk number 20: Multivariate_Extremes.Rnw:233-237
###################################################
m1 <- fbvpot(loss, thresh, model = "alog", asy1 = 1)
m2 <- fbvpot(loss, thresh, model = "bilog")
m3 <- fbvpot(loss, thresh, model = "bilog", likelihood = "poisson")
round(rbind(fitted(m2), std.errors(m2)), 3)


###################################################
### code chunk number 21: Multivariate_Extremes.Rnw:242-247
###################################################
abvnonpar(data = loss, method = "pot", k = k0, epmar = TRUE, 
  plot = TRUE, lty = 3)
plot(m1, which = 2, add = TRUE)
plot(m2, which = 2, add = TRUE, lty = 4)
plot(m3, which = 2, add = TRUE, lty = 2)


###################################################
### code chunk number 22: qcthresh
###################################################
lts <- c(1e-04, 100)
plot(loss, log = "xy", col = "grey", xlim = lts, ylim = lts)
plot(m1, which = 3, p = c(0.95,0.975,0.99), tlty = 0, add = TRUE)
abline(v=thresh[1], h=thresh[2])


###################################################
### code chunk number 23: Multivariate_Extremes.Rnw:261-262
###################################################
lts <- c(1e-04, 100)
plot(loss, log = "xy", col = "grey", xlim = lts, ylim = lts)
plot(m1, which = 3, p = c(0.95,0.975,0.99), tlty = 0, add = TRUE)
abline(v=thresh[1], h=thresh[2])


###################################################
### code chunk number 24: chiplot
###################################################
old <- par(mfrow = c(2,1))
chiplot(loss, ylim1 = c(-0.25,1), ylim2 = c(-0.25,1), nq = 200, 
  qlim = c(0.02,0.98), which = 1:2, spcases = TRUE)
par(old)


###################################################
### code chunk number 25: Multivariate_Extremes.Rnw:285-286
###################################################
old <- par(mfrow = c(2,1))
chiplot(loss, ylim1 = c(-0.25,1), ylim2 = c(-0.25,1), nq = 200, 
  qlim = c(0.02,0.98), which = 1:2, spcases = TRUE)
par(old)


###################################################
### code chunk number 26: etaplot
###################################################
fla <- apply(-1/log(ula), 1, min)
thresh <- quantile(fla, probs = c(0.025, 0.975))
tcplot(fla, thresh, nt = 100, pscale = TRUE, which = 2, vci = FALSE, 
  cilty = 2, type = "l", ylim = c(-0.2,1.2), ylab = "Tail Dependence")
abline(h = c(0,1))


###################################################
### code chunk number 27: Multivariate_Extremes.Rnw:304-307
###################################################
thresh <- quantile(fla, probs = 0.8)
m1 <- fpot(fla, thresh = thresh)
cat("Tail Dependence:", fitted(m1)["shape"], "\n")


###################################################
### code chunk number 28: Multivariate_Extremes.Rnw:310-312
###################################################
m2 <- fpot(fla, thresh = thresh, shape = 1)
anova(m1, m2, half = TRUE)


###################################################
### code chunk number 29: Multivariate_Extremes.Rnw:317-318
###################################################
fla <- apply(-1/log(ula), 1, min)
thresh <- quantile(fla, probs = c(0.025, 0.975))
tcplot(fla, thresh, nt = 100, pscale = TRUE, which = 2, vci = FALSE, 
  cilty = 2, type = "l", ylim = c(-0.2,1.2), ylab = "Tail Dependence")
abline(h = c(0,1))