1"geweke.diag" <-  function (x, frac1 = 0.1, frac2 = 0.5)
2{
3    if (frac1 < 0 || frac1 > 1) {
4        stop("frac1 invalid")
5    }
6    if (frac2 < 0 || frac2 > 1) {
7        stop("frac2 invalid")
8    }
9    if (frac1 + frac2 > 1) {
10        stop("start and end sequences are overlapping")
11    }
12    if (is.mcmc.list(x)) {
13        return(lapply(x, geweke.diag, frac1, frac2))
14    }
15    x <- as.mcmc(x)
16    xstart <- c(start(x), floor(end(x) - frac2 * (end(x) - start(x))))
17    xend <- c(ceiling(start(x) + frac1 * (end(x) - start(x))), end(x))
18    y.variance <- y.mean <- vector("list", 2)
19    for (i in 1:2) {
20        y <- window(x, start = xstart[i], end = xend[i])
21        y.mean[[i]] <- apply(as.matrix(y), 2, mean)
22        y.variance[[i]] <- spectrum0.ar(y)$spec/niter(y)
23    }
24    z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
25    out <- list(z = z, frac = c(frac1, frac2))
26    class(out) <- "geweke.diag"
27    return(out)
28}
29
30"geweke.plot" <-
31  function (x, frac1 = 0.1, frac2 = 0.5, nbins = 20,
32            pvalue = 0.05, auto.layout = TRUE, ask, ...)
33{
34  if (missing(ask)) {
35    ask <- if (is.R()) {
36      dev.interactive()
37    }
38    else {
39      interactive()
40    }
41  }
42
43  x <- as.mcmc.list(x)
44  oldpar <- NULL
45  on.exit(par(oldpar))
46  if (auto.layout)
47    oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),
48                    Nparms = nvar(x)))
49  ystart <- seq(from = start(x), to = (start(x) + end(x))/2, length = nbins)
50  if (is.R())
51    gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)),
52               dimnames = list(ystart, varnames(x), chanames(x)))
53  else
54    gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)),
55               dimnames = list(ystart, varnames(x), chanames(x)))
56
57  for (n in 1:length(ystart)) {
58    geweke.out <- geweke.diag(window(x, start = ystart[n]),
59                              frac1 = frac1, frac2 = frac2)
60    for (k in 1:nchain(x)) gcd[n, , k] <- geweke.out[[k]]$z
61  }
62  climit <- qnorm(1 - pvalue/2)
63  for (k in 1:nchain(x)) for (j in 1:nvar(x)) {
64    ylimit <- max(c(climit, abs(gcd[, j, k])))
65    plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment",
66         ylab = "Z-score", pch = 4, ylim = c(-ylimit, ylimit),
67         ...)
68    abline(h = c(climit, -climit), lty = 2)
69    if (nchain(x) > 1) {
70      title(main = paste(varnames(x, allow.null = FALSE)[j],
71              " (", chanames(x, allow.null = FALSE)[k], ")",
72              sep = ""))
73    }
74    else {
75      title(main = paste(varnames(x, allow.null = FALSE)[j],
76              sep = ""))
77    }
78    if (k==1 && j==1)
79       oldpar <- c(oldpar, par(ask = ask))
80  }
81  invisible(list(start.iter = ystart, z = gcd))
82}
83
84"print.geweke.diag" <- function (x, digits = min(4, .Options$digits), ...)
85  ## Print method for output from geweke.diag
86{
87  cat("\nFraction in 1st window =", x$frac[1])
88  cat("\nFraction in 2nd window =", x$frac[2], "\n\n")
89  print.default(x$z, digits = digits, ...)
90  cat("\n")
91  invisible(x)
92}
93
94
95
96
97
98
99
100
101
102
103
104