normalCopula(c(0.5, 0.6, 0.7), dim = 3, dispstr = "un")
t.cop <- tCopula(c(0.5, 0.3), dim = 3, dispstr = "toep",
df = 2, df.fixed = TRUE)
getSigma(t.cop) # P matrix (with diagonal = 1)
stopifnot(all.equal(toeplitz(c(1, .5, .3)), getSigma(t.cop)))
## dispersion "AR1" :
nC.7 <- normalCopula(0.8, dim = 7, dispstr = "ar1")
getSigma(nC.7)
stopifnot(all.equal(toeplitz(.8^(0:6)), getSigma(nC.7)))
## from the wrapper
norm.cop <- ellipCopula("normal", param = c(0.5, 0.6, 0.7),
dim = 3, dispstr = "un")
if(require("scatterplot3d") && dev.interactive(orNone=TRUE)) {
## 3d scatter plot of 1000 random observations
scatterplot3d(rCopula(1000, norm.cop))
scatterplot3d(rCopula(1000, t.cop))
}
set.seed(12)
uN <- rCopula(512, norm.cop)
set.seed(2); pN1 <- pCopula(uN, norm.cop)
set.seed(3); pN2 <- pCopula(uN, norm.cop)
stopifnot(identical(pN1, pN2)) # no longer random for dim = 3
(Xtras <- copula:::doExtras())
if(Xtras) { ## a bit more accurately:
set.seed(4); pN1. <- pCopula(uN, norm.cop, abseps = 1e-9)
set.seed(5); pN2. <- pCopula(uN, norm.cop, abseps = 1e-9)
stopifnot(all.equal(pN1., pN2., 1e-5))# see 3.397e-6
## but increasing the required precision (e.g., abseps=1e-15) does *NOT* help
}
## For smaller copula dimension 'd', alternatives are available and
## non-random, see ?GenzBretz from package 'mvtnorm' :
has_mvtn <- "package:mvtnorm" %in% search() #% << (workaround ESS Rd render bug)
if(!has_mvtn)
require("mvtnorm")# -> GenzBretz(), Miva(), and TVPACK() are available
## Note that Miwa() would become very slow for dimensions 5, 6, ..
set.seed(4); pN1.M <- pCopula(uN, norm.cop, algorithm = Miwa(steps = 512))
set.seed(5); pN2.M <- pCopula(uN, norm.cop, algorithm = Miwa(steps = 512))
stopifnot(all.equal(pN1.M, pN2.M, tol= 1e-15))# *no* randomness
set.seed(4); pN1.T <- pCopula(uN, norm.cop, algorithm = TVPACK(abseps = 1e-10))
set.seed(5); pN2.T <- pCopula(uN, norm.cop, algorithm = TVPACK(abseps = 1e-14))
stopifnot(all.equal(pN1.T, pN2.T, tol= 1e-15))# *no* randomness (but no effect of 'abseps')
if(!has_mvtn)
detach("package:mvtnorm")# (revert)
## Versions with unspecified parameters:
tCopula()
allEQ <- function(u,v) all.equal(u, v, tolerance=0)
stopifnot(allEQ(ellipCopula("norm"), normalCopula()),
allEQ(ellipCopula("t"), tCopula()))
tCopula(dim=3)
tCopula(dim=4, df.fixed=TRUE)
tCopula(dim=5, disp = "toep", df.fixed=TRUE)
normalCopula(dim=4, disp = "un")
## Toeplitz after *permutation* dispersions (new in copula 1.1-0) ---------
tpar <- c(7,5,3)/8 # *gives* pos.def.:
(ev <- eigen(toeplitz(c(1, tpar)), symmetric=TRUE, only.values=TRUE)$values)
stopifnot(ev > 0)
N4. <- ellipCopula("normal", dim=4, param=tpar, dispstr = "toep") #"regular"
## reversed order is "the same" for toeplitz structure:
N4.pr <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(4:1))
N4.p1 <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(c(4,1:3)))
N4.p2 <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(c(4:3,1:2)))
N4.p3 <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(c(2,4,1,3)))
(pm <- attr(N4.p3@dispstr, "perm")) # (2 4 1 3)
ip <- c(3,1,4,2) # the *inverse* permutation of (2 4 1 3) = Matrix::invPerm(pm)
(Sp3 <- getSigma(N4.p3)) # <-- "permuted toeplitz"
Sp3[ip, ip] # re-ordered rows & columns => *is* toeplitz :
stopifnot(exprs = {
## permutations pm and ip are inverses:
pm[ip] == 1:4
ip[pm] == 1:4
is.matrix(T4 <- toeplitz(c(1, tpar)))
identical(getSigma(N4.), T4)
identical(getSigma(N4.pr), T4) # 4:1 and 1:4 is "the same" for Rho
identical(Sp3[ip, ip] , T4)
identical(Sp3, T4[pm,pm])
})
## Data generation -- NB: The U matrices are equal only "in distribution":
set.seed(7); U.p3 <- rCopula(1000, N4.p3)
set.seed(7); U. <- rCopula(1000, N4.)
stopifnot(exprs = {
all.equal(loglikCopula(tpar, u=U.p3, copula= N4.p3),
loglikCopula(tpar, u=U.p3[,ip], copula= N4.) -> LL3)
all.equal(loglikCopula(tpar, u=U., copula= N4.),
loglikCopula(tpar, u=U.[,pm], copula= N4.p3) -> LL.)
})
c(LL. , LL3)# similar but different
if(Xtras) {
fm. <- fitCopula(N4. , U. )
fm.3 <- fitCopula(N4.p3, U.p3)
summary(fm.3)
stopifnot(all.equal(coef(fm.), coef(fm.3), tol = 0.01))# similar but different
}
Run the code above in your browser using DataLab