mcmc (version 0.7-2)

initseq: Initial Sequence Estimators

Description

Variance of sample mean of functional of reversible Markov chain using methods of Geyer (1992).

Usage

initseq(x)

Arguments

x

a numeric vector that is a scalar-valued functional of a reversible Markov chain.

Value

• a list containing the following components:
• gamma0
  the scalar $\gamma_0$, the marginal variance of x.
• Gamma.pos
  the vector $\Gamma$, estimated so as to be nonnegative, where, as always, Ruses one-origin indexing so Gamma.pos[1] is $\Gamma_0$.
• Gamma.dec
  the vector $\Gamma$, estimated so as to be nonnegative and nonincreasing, where, as always, Ruses one-origin indexing so Gamma.dec[1] is $\Gamma_0$.
• Gamma.con
  the vector $\Gamma$, estimated so as to be nonnegative and nonincreasing and convex, where, as always, Ruses one-origin indexing so Gamma.con[1] is $\Gamma_0$.
• var.pos
  the scalar - gamma0 + 2 * sum(Gamma.pos), which is the asymptotic variance in the Markov chain CLT. Divide by length(x) to get the approximate variance of the sample mean of x.
• var.dec
  the scalar - gamma0 + 2 * sum(Gamma.dec), which is the asymptotic variance in the Markov chain CLT. Divide by length(x) to get the approximate variance of the sample mean of x.
• var.con
  the scalar - gamma0 + 2 * sum(Gamma.con), which is the asymptotic variance in the Markov chain CLT. Divide by length(x) to get the approximate variance of the sample mean of x.

Details

Let $$\gamma_k = \textrm{cov}(X_i, X_{i + k})$$ considered as a function of the lag $k$ be the autocovariance function of the input time series. Define $$\Gamma_k = \gamma_{2 k} + \gamma_{2 k + 1}$$ the sum of consecutive pairs of autocovariances. Then Theorem 3.1 in Geyer (1992) says that $\Gamma_k$ considered as a function of $k$ is strictly positive, strictly decreasing, and strictly convex, assuming the input time series is a scalar-valued functional of a reversible Markov chain. All of the MCMC done by this package is reversible. This Rfunction estimates the big gamma function, $\Gamma_k$ considered as a function of $k$, subject to three different constraints, (1) nonnegative, (2) nonnegative and nonincreasing, and (3) nonnegative, nonincreasing, and convex. It also estimates the variance in the Markov chain central limit theorem (CLT) $$\gamma_0 + 2 \sum_{k = 1}^\infty \gamma_k = - \gamma_0 + 2 \sum_{k = 0}^\infty \Gamma_k$$

Note: The batch means provided by metrop are also scalar functionals of a reversible Markov chain. Thus these initial sequence estimators applied to the batch means give valid standard errors for the mean of the match means even when the batch length is too short to provide a valid estimate of asymptotic variance. One does, of course, have to multiply the asymptotic variance of the batch means by the batch length to get the asymptotic variance for the unbached chain.

References

Geyer, C. J. (1992) Practical Markov Chain Monte Carlo. Statistical Science 7 473--483.

See Also

metrop

Examples

n <- 2e4
rho <- 0.99
x <- arima.sim(model = list(ar = rho), n = n)
out <- initseq(x)
plot(seq(along = out$Gamma.pos) - 1, out$Gamma.pos,
   xlab = "k", ylab = expression(Gamma[k]), type = "l")
lines(seq(along = out$Gamma.dec) - 1, out$Gamma.dec, col = "red")
lines(seq(along = out$Gamma.con) - 1, out$Gamma.con, col = "blue")
# asymptotic 95\% confidence interval for mean of x
mean(x) + c(-1, 1) * qnorm(0.975) * sqrt(out$var.con / length(x))
# estimated asymptotic variance
out$var.con
# theoretical asymptotic variance
(1 + rho) / (1 - rho) * 1 / (1 - rho^2)
# illustrating use with batch means
bm <- apply(matrix(x, nrow = 5), 2, mean)
initseq(bm)$var.con * 5