varianceReduction | R Documentation |
Computing antithetic variates or Latin hypercube samples.
rAntitheticVariates(u) rLatinHypercube(u, ...)
u |
a n x d- |
... |
additional arguments passed to the underlying |
rAntitheticVariates()
takes any copula sample u, builds
1-u, and returns the two matrices in the form of an array; this
can be used for the variance-reduction method of (componentwise)
antithetic variates.
rLatinHypercube()
takes any copula sample,
shifts its marginal ranks minus 1 by standard uniforms and then
divides by the sample size in order to obtain a Latin hypercubed
sample.
rAntitheticVariates()
array
of dimension
n x d x 2, say r
, where r[,,1]
contains the original sample u
and r[,,2]
contains
the sample 1-u
.
rLatinHypercube()
matrix
of the same
dimensions as u
.
Cambou, M., Hofert, M. and Lemieux, C. (2016). Quasi-random numbers for copula models. Statistics and Computing, 1–23.
Packham, N. and Schmidt, W. M. (2010). Latin hypercube sampling with dependence and applications in finance. Journal of Computational Finance 13(3), 81–111.
### 1 Basic plots ############################################################## ## Generate data from a Gumbel copula cop <- gumbelCopula(iTau(gumbelCopula(), tau = 0.5)) n <- 1e4 set.seed(271) U <- rCopula(n, copula = cop) ## Transform the sample to a Latin Hypercube sample U.LH <- rLatinHypercube(U) ## Plot ## Note: The 'variance-reducing property' is barely visible, but that's okay layout(rbind(1:2)) plot(U, xlab = quote(U[1]), ylab = quote(U[2]), pch = ".", main = "U") plot(U.LH, xlab = quote(U[1]), ylab = quote(U[2]), pch = ".", main = "U.LH") layout(1) # reset layout ## Transform the sample to an Antithetic variate sample U.AV <- rAntitheticVariates(U) stopifnot(identical(U.AV[,,1], U ), identical(U.AV[,,2], 1-U)) ## Plot original sample and its corresponding (componentwise) antithetic variates layout(rbind(1:2)) plot(U.AV[,,1], xlab = quote(U[1]), ylab = quote(U[2]), pch=".", main= "U") plot(U.AV[,,2], xlab = quote(U[1]), ylab = quote(U[2]), pch=".", main= "1 - U") layout(1) # reset layout ### 2 Small variance-reduction study for exceedance probabilities ############## ## Auxiliary function for approximately computing P(U_1 > u_1,..., U_d > u_d) ## by Monte Carlo simulation based on pseudo-random numbers, Latin hypercube ## sampling and quasi-random numbers. sProb <- function(n, copula, u) { d <- length(u) stopifnot(n >= 1, inherits(copula, "Copula"), 0 < u, u < 1, d == dim(copula)) umat <- rep(u, each = n) ## Pseudo-random numbers U <- rCopula(n, copula = copula) PRNG <- mean(rowSums(U > umat) == d) ## Latin hypercube sampling (based on the recycled 'U') U. <- rLatinHypercube(U) LHS <- mean(rowSums(U. > umat) == d) ## Quasi-random numbers U.. <- cCopula(sobol(n, d = d, randomize = TRUE), copula = copula, inverse = TRUE) QRNG <- mean(rowSums(U.. > umat) == d) ## Return c(PRNG = PRNG, LHS = LHS, QRNG = QRNG) } ## Simulate the probabilities of falling in (u_1,1] x ... x (u_d,1] library(qrng) # for quasi-random numbers (Xtras <- copula:::doExtras()) # determine whether examples will be extra (long) B <- if(Xtras) 500 else 100 # number of replications n <- if(Xtras) 1000 else 200 # sample size d <- 2 # dimension; note: for d > 2, the true value depends on the seed nu <- 3 # degrees of freedom th <- iTau(tCopula(df = nu), tau = 0.5) # correlation parameter cop <- tCopula(param = th, dim = d, df = nu) # t copula u <- rep(0.99, d) # lower-left endpoint of the considered cube set.seed(42) # for reproducibility true <- prob(cop, l = u, u = rep(1, d)) # true exceedance probability system.time(res <- replicate(B, sProb(n, copula = cop, u = u))) ## "abbreviations": PRNG <- res["PRNG",] LHS <- res["LHS" ,] QRNG <- res["QRNG",] ## Compute the variance-reduction factors and % improvements vrf <- var(PRNG) / var(LHS) # variance reduction factor w.r.t. LHS vrf. <- var(PRNG) / var(QRNG) # variance reduction factor w.r.t. QRNG pim <- (var(PRNG) - var(LHS)) / var(PRNG) *100 # improvement w.r.t. LHS pim. <- (var(PRNG) - var(QRNG))/ var(PRNG) *100 # improvement w.r.t. QRNG ## Boxplot boxplot(list(PRNG = PRNG, LHS = LHS, QRNG = QRNG), notch = TRUE, main = substitute("Simulated exceedance probabilities" ~ P(bold(U) > bold(u))~~ "for a" ~ t[nu.]~"copula", list(nu. = nu)), sub = sprintf( "Variance-reduction factors and %% improvements: %.1f (%.0f%%), %.1f (%.0f%%)", vrf, pim, vrf., pim.)) abline(h = true, lty = 3) # true value mtext(sprintf("B = %d replications with n = %d and d = %d", B, n, d), side = 3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.