Nothing
BiCopGofKendall <- function(u1, u2, family, B = 100, level = 0.05) {
if (is.null(u1) == TRUE || is.null(u2) == TRUE)
stop("u1 and/or u2 are not set or have length zero.")
if (any(u1 > 1) || any(u1 < 0))
stop("Data has be in the interval [0,1].")
if (any(u2 > 1) || any(u2 < 0))
stop("Data has be in the interval [0,1].")
if (length(u1) != length(u2))
stop("Lengths of 'u1' and 'u2' do not match.")
if (length(u1) < 2)
stop("Number of observations has to be at least 2.")
if (!(family %in% c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 13, 14, 16, 17, 18, 19, 20, 23, 24, 26, 27, 28,
29, 30, 33, 34, 36, 37, 38, 39, 40)))
stop("Copula family not implemented.")
if (level < 0 & level > 1)
stop("Significance level has to be between 0 and 1.")
if (family %in% c(13, 14, 16, 17, 18, 19, 20)) {
u1 <- 1 - u1
u2 <- 1 - u2
family <- family - 10
} else if (family %in% c(23, 24, 26, 27, 28, 29, 30)) {
u1 <- 1 - u1
family <- family - 20
} else if (family %in% c(33, 34, 36, 37, 38, 39, 40)) {
u2 <- 1 - u2
family <- family - 30
}
ostat <- obs.stat(u1, u2, family)
if (B == 0) {
# no bootstrap
sn.obs <- ostat$Sn
tn.obs <- ostat$Tn
out <- list(Sn = sn.obs, Tn = tn.obs)
} else {
bstat <- list()
for (i in 1:B) bstat[[i]] <- boot.stat(u1, u2, family)
sn.boot <- rep(0, B)
tn.boot <- rep(0, B)
for (i in 1:B) {
sn.boot[i] <- bstat[[i]]$sn
tn.boot[i] <- bstat[[i]]$tn
}
sn.obs <- ostat$Sn
tn.obs <- ostat$Tn
k <- as.integer((1 - level) * B)
sn.critical <- sn.boot[k] # critical value of test at level 0.05
tn.critical <- tn.boot[k] # critical value of test at level 0.05
pv.sn <- sapply(sn.obs, function(x) (1/B) * length(which(sn.boot[1:B] >= x))) # P-value of Sn
pv.tn <- sapply(tn.obs, function(x) (1/B) * length(which(tn.boot[1:B] >= x))) # P-value of Tn
out <- list(p.value.CvM = pv.sn, p.value.KS = pv.tn, statistic.CvM = sn.obs, statistic.KS = tn.obs)
}
return(out)
}
boot.stat <- function(u, v, fam) {
n <- length(u)
t <- seq(1, n)/(n + 1e-04)
kt <- rep(0, n)
# estimate paramemter for different copula family from (u,v)
param <- suppressWarnings({
BiCopEst(u, v, family = fam)
})
# calulate k(t) and kn(t) of bootstrap sample data
if (fam == 1) {
# normal
sam <- BiCopSim(n, 1, param$par, param$par2) # generate data for the simulation of K(t)
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam) # parameter estimation of sample data
sim <- BiCopSim(10000, 1, sam.par$par, sam.par$par2) # generate data for the simulation of theo. K(t)
cormat <- matrix(c(1, param$par, param$par, 1), 2, 2)
dcop <- rep(0, 10000)
for (i in 1:10000) dcop[i] <- pmvnorm(upper = c(qnorm(sim[i, 1]), qnorm(sim[i, 2])), corr = cormat)
kt <- sapply(t, function(x) (1/10000) * length(which(dcop[1:10000] <= x))) # simulate K(t) of sample data
}
if (fam == 2) {
# t
sam <- BiCopSim(n, fam, param$par, param$par2) # generate data for the simulation of K(t)
sam.par <- suppressWarnings({
BiCopEst(sam[, 1], sam[, 2], family = fam)
}) # parameter estimation of sample data
sim <- BiCopSim(10000, fam, sam.par$par, sam.par$par2) # generate data for the simulation of theo. K(t)
cormat <- matrix(c(1, param$par, param$par, 1), 2, 2)
dcop <- rep(0, 10000)
for (i in 1:10000) dcop[i] <- pmvt(upper = c(qt(sim[i, 1], df = param$par2), qt(sim[i, 2], df = param$par2)),
corr = cormat, df = param$par2)
kt <- sapply(t, function(x) (1/10000) * length(which(dcop[1:10000] <= x))) # simulate K(t) of sample data
} else if (fam == 3) {
# Clayton
sam <- BiCopSim(n, 3, param$par) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam)$par # estimate parameter of sample data
kt <- t + t * (1 - t^sam.par)/sam.par
} else if (fam == 4) {
# gumbel
sam <- BiCopSim(n, 4, param$par) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam)$par # estimate parameter of sample data
kt <- t - t * log(t)/(sam.par)
} else if (fam == 5) {
# frank
sam <- BiCopSim(n, 5, param$par) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam)$par # estimate parameter of sample data
kt <- t + log((1 - exp(-sam.par))/(1 - exp(-sam.par * t))) * (1 - exp(-sam.par * t))/(sam.par * exp(-sam.par *
t))
} else if (fam == 6) {
# joe
sam <- BiCopSim(n, 6, param$par) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam)$par # estimate parameter of sample data
kt <- t - (log(1 - (1 - t)^sam.par) * (1 - (1 - t))^sam.par)/(sam.par * (1 - t)^(sam.par - 1))
} else if (fam == 7) {
# BB1
sam <- BiCopSim(n, 7, param$par, param$par2) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam) # estimate parameter of sample data
theta <- sam.par$par
delta <- sam.par$par2
kt <- t + 1/(theta * delta) * (t^(-theta) - 1)/(t^(-1 - theta))
} else if (fam == 8) {
# BB6
sam <- BiCopSim(n, 8, param$par, param$par2) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam) # estimate parameter of sample data
theta <- sam.par$par
delta <- sam.par$par2
kt <- t + log(-(1 - t)^theta + 1) * (1 - t - (1 - t)^(-theta) + (1 - t)^(-theta) * t)/(delta * theta)
} else if (fam == 9) {
# BB7
sam <- BiCopSim(n, 9, param$par, param$par2) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam) # estimate parameter of sample data
theta <- sam.par$par
delta <- sam.par$par2
kt <- t + 1/(theta * delta) * ((1 - (1 - t)^theta)^(-delta) - 1)/((1 - t)^(theta - 1) * (1 - (1 -
t)^theta)^(-delta - 1))
} else if (fam == 10) {
# BB8
sam <- BiCopSim(n, 10, param$par, param$par2) # generate sample data
sam.par <- BiCopEst(sam[, 1], sam[, 2], family = fam) # estimate parameter of sample data
theta <- sam.par$par
delta <- sam.par$par2
kt <- t + log(((1 - t * delta)^theta - 1)/((1 - delta)^theta - 1)) * (1 - t * delta - (1 - t * delta)^(-theta) +
(1 - t * delta)^(-theta) * t * delta)/(theta * delta)
}
# calculate emp. Kn
w <- rep(0, n)
w[1:n] <- mapply(function(x, y) (1/n) * length(which(x > sam[, 1] & y > sam[, 2])), sam[, 1], sam[, 2])
w <- sort(w)
kn <- rep(0, n)
kn <- sapply(t, function(x) (1/n) * length(which(w[1:n] <= x)))
# calculate test statistic Sn
Sn1 <- 0
Sn2 <- 0
for (j in 1:(n - 1)) {
Sn1 <- Sn1 + ((kn[j])^2 * (kt[j + 1] - kt[j]))
Sn2 <- Sn2 + (kn[j]) * ((kt[j + 1])^2 - (kt[j])^2)
}
sn <- n/3 + n * Sn1 - n * Sn2
# calculation of test statistics Tn
tm <- matrix(0, n - 1, 2)
# mit i=0
for (j in 1:(n - 1)) {
tm[j, 1] <- abs(kn[j] - kt[j])
}
# mit i=1
for (j in 1:(n - 1)) {
tm[j, 2] <- abs(kn[j] - kt[j + 1])
}
tn <- max(tm) * sqrt(n)
sn <- sort(sn) # vector of ordered statistic Sn
tn <- sort(tn) # vector of ordered statistic Tn
out <- list(sn = sn, tn = tn)
}
obs.stat <- function(u, v, fam) {
n <- length(u)
t <- seq(1, n)/(n + 1e-04)
kt <- rep(0, n)
# estimate paramemter for different copula family from (u,v)
param <- suppressWarnings({
BiCopEst(u, v, family = fam)
})
# calculate observed K(t) of (u,v)
kt.obs <- rep(0, n)
if (fam == 1) {
sim <- BiCopSim(10000, 1, param$par) # generate data for the simulation of K(t)
cormat <- matrix(c(1, param$par, param$par, 1), 2, 2)
dcop <- rep(0, 10000)
for (i in 1:10000) dcop[i] <- pmvnorm(upper = c(qnorm(sim[i, 1]), qnorm(sim[i, 2])), corr = cormat)
kt.obs <- sapply(t, function(x) (1/10000) * length(which(dcop[1:10000] <= x))) # simulate K(t) of sample data
} else if (fam == 2) {
sim <- BiCopSim(10000, 2, param$par, param$par2) # generate data for the simulation of K(t)
cormat <- matrix(c(1, param$par, param$par, 1), 2, 2)
dcop <- rep(0, 10000)
for (i in 1:10000) dcop[i] <- pmvt(upper = c(qt(sim[i, 1], df = param$par2), qt(sim[i, 2], df = param$par2)),
corr = cormat, df = param$par2)
kt.obs <- sapply(t, function(x) (1/10000) * length(which(dcop[1:10000] <= x))) # simulate K(t) of sample data
} else if (fam == 3) {
kt.obs <- t + t * (1 - t^param$par)/param$par
} else if (fam == 4) {
kt.obs <- t - t * log(t)/(param$par)
} else if (fam == 5) {
kt.obs <- t + log((1 - exp(-param$par))/(1 - exp(-param$par * t))) * (1 - exp(-param$par * t))/(param$par *
exp(-param$par * t))
} else if (fam == 6) {
kt.obs <- t - (log(1 - (1 - t)^param$par) * (1 - (1 - t))^param$par)/(param$par * (1 - t)^(param$par -
1))
} else if (fam == 7) {
theta <- param$par
delta <- param$par2
kt.obs <- t + 1/(theta * delta) * (t^(-theta) - 1)/(t^(-1 - theta))
} else if (fam == 8) {
theta <- param$par
delta <- param$par2
kt.obs <- t + log(-(1 - t)^theta + 1) * (1 - t - (1 - t)^(-theta) + (1 - t)^(-theta) * t)/(delta *
theta)
} else if (fam == 9) {
theta <- param$par
delta <- param$par2
kt.obs <- t + 1/(theta * delta) * ((1 - (1 - t)^theta)^(-delta) - 1)/((1 - t)^(theta - 1) * (1 -
(1 - t)^theta)^(-delta - 1))
} else if (fam == 10) {
theta <- param$par
delta <- param$par2
kt.obs <- t + log(((1 - t * delta)^theta - 1)/((1 - delta)^theta - 1)) * (1 - t * delta - (1 - t *
delta)^(-theta) + (1 - t * delta)^(-theta) * t * delta)/(theta * delta)
}
# calculation of observed Kn
w <- rep(0, n)
w[1:n] <- mapply(function(x, y) (1/n) * length(which(x > u & y > v)), u, v)
w <- sort(w)
kn.obs <- rep(0, n)
kn.obs <- sapply(t, function(x) (1/n) * length(which(w[1:n] <= x)))
# calculation of observed value Sn
Sn1 <- 0
Sn2 <- 0
for (j in 1:(n - 1)) {
Sn1 <- Sn1 + ((kn.obs[j])^2 * (kt.obs[j + 1] - kt.obs[j]))
Sn2 <- Sn2 + (kn.obs[j]) * ((kt.obs[j + 1])^2 - (kt.obs[j])^2)
}
Sn <- n/3 + n * Sn1 - n * Sn2 # observed Sn
# calculation of observed Tn
tn.obs <- matrix(0, n - 1, 2)
# mit i=0
for (j in 1:(n - 1)) {
tn.obs[j, 1] <- abs(kn.obs[j] - kt.obs[j])
}
# mit i=1
for (j in 1:(n - 1)) {
tn.obs[j, 2] <- abs(kn.obs[j] - kt.obs[j + 1])
}
Tn <- max(tn.obs) * sqrt(n)
out <- list(Sn = Sn, Tn = Tn)
return(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.