tau <- 0.5
(theta <- copGumbel@iTau(tau)) # 2
d <- 20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))
## basic check empirical Kendall distribution function
set.seed(271)
n <- 1000
U <- rCopula(n, cop)
X <- qnorm(U)
K.sample <- pCopula(U, copula=cop)
u <- seq(0, 1, length.out=256)
edfK <- ecdf(K.sample)
plot(u, edfK(u), type="l", ylim=c(0,1),
xlab=expression(italic(u)), ylab=expression(K[n](italic(u)))) # simulated
K.n <- Kn(u, x=X)
lines(u, K.n, col="royalblue3") # Kn
## difference at 0
edfK(0) # edf of K at 0
K.n[1] # K_n(0); this is > 0 since K.n is the edf of a discrete distribution
## => therefore, Kn(K.sample, x=X) is not uniform
plot(Kn(K.sample, x=X), ylim=c(0,1))
## Note: Kn(0) -> 0 for n -> Inf
## compute Kendall distribution function
u <- seq(0,1, length.out = 255)
Ku <- pK(u, cop=cop@copula, d=d) # exact
Ku.MC <- pK(u, cop=cop@copula, d=d, n.MC=1000) # via Monte Carlo
## build sample from K
set.seed(1)
n <- 200
W <- rK(n, cop)
## plot empirical distribution function based on W
## and the corresponding theoretical Kendall distribution function
## (exact and via Monte Carlo)
plot(ecdf(W), col="blue", xlim=c(0,1), verticals=TRUE,
main = expression("Empirical"~ F[n]( C(U) ) ~
"and its Kendall distribution"~ K(u)),
do.points=FALSE, asp=1)
abline(0,1, lty=2); abline(h=0:1, v=0:1, lty=3, col="gray")
lines(u, Ku.MC, col="red")
lines(u, Ku, col="black")
legend(.2,.6, expression(F[n],K(u), K[MC](u)),
col=c("blue","red","black"), lty=1, bty="n",
xjust = 0.25, yjust = 0.5)
### testing qK
uexpr <- quote( 0:63/63 ); u <- eval(uexpr)
d <- 10
cop <- onacopulaL("Gumbel", list(theta = 2, 1:d))
system.time(Ku1 <- qK(u, cop=cop@copula, d=d, method="simple"))
system.time(Ku2 <- qK(u, cop=cop@copula, d=d, method="sort"))
system.time(Ku3 <- qK(u, cop=cop@copula, d=d, method="discrete",
u.grid=0:1e4/1e4))
system.time(Ku4 <- qK(u, cop=cop@copula, d=d, method="monoH.FC",
u.grid=0:5e2/5e2))
cols <- adjustcolor(c("gray80", "light blue", "royal blue", "purple3"), 0.6)
matplot(u, cbind(Ku1,Ku2,Ku3,Ku4), type="l", lwd=2*4:1, lty = 1:4, col= cols,
xlab=substitute(u == U, list(U=uexpr)), ylab=expression({K^{-1}}(u)))
legend("topleft", col=cols, lwd=2*4:1, lty=1:4, bty="n", inset=.03,
legend= paste0("method= ",
sQuote(c("simple", "sort", "discrete", "monoH.FC"))))
Run the code above in your browser using DataLab