R: Multivariate Monte Carlo Standard Errors Estimation.
mcse.multi
R Documentation
Multivariate Monte Carlo Standard Errors Estimation.
Description
Function returns the estimate of the covariance matrix in the Markov Chain CLT using batch means or spectral variance methods (with different lag windows). The function also returns the volume of the resulting ellipsoidal confidence regions.
Usage
mcse.multi(x, method = "bm", size = "sqroot", g = NULL, level = 0.95, large = FALSE)
Arguments
x
an n by p matrix that of the Markov chain output
method
any of ``bm'', ``bartlett'', ``tukey''. ``bm'' represents batch means estimator, ``bartlett'' and ``tukey'' represents the modified-Bartlett window and the Tukey-Hanning windows for the spectral variance estimators.
size
can take character values of ``sqroot'' and ``cuberoot'' or any numeric value between 1 and n. Size represents the batch size in bm and the truncation point in bartlett and tukey. sqroot means size is floor(n^(1/2) and cuberoot means size is floor(n^(1/3)).
g
a function that represents features of interest. g is applied to each row of x and thus g should take a vector input only. If g is NULL, g is set to be identity, which is estimation of the mean of the target density.
level
confidence level of the confidence ellipsoid.
large
if TRUE, returns the volume of the large sample confidence region using a chi square critical value.
Value
A list is returned with the following components,
cov
a p by p covariance matrix estimate.
vol
volume of the confidence ellipsoid to the pth root.
est
estimate of g(x).
nsim
number of rows of the input x.
method
method used to calculate matrix cov.
large
logical of whether a large sample confidence region volume was calculated.
size
value of size used to calculate cov.
See Also
mcse, which acts on a vector.
mcse.mat, which applies mcse to each
column of a matrix or data frame.
mcse.q and mcse.q.mat, which
compute standard errors for quantiles.
Examples
library(mAr)
p <- 3
n <- 1e3
omega <- 5*diag(1,p)
## Making correlation matrix var(1) model
set.seed(100)
foo <- matrix(rnorm(p^2), nrow = p)
foo <- foo %*% t(foo)
phi <- foo / (max(eigen(foo)$values) + 1)
out <- as.matrix(mAr.sim(rep(0,p), phi, omega, N = n))
mcse.bm <- mcse.multi(x = out)
mcse.tuk <- mcse.multi(x = out, method = "tukey")
# If we are only estimating the mean of the first component,
# and the second moment of the second component
g <- function(x) return(c(x[1], x[2]^2))
mcse <- mcse.multi(x = out, g = g)
providing 
.