1
2############## Pesaran's CD test and Breusch/Pagan LM Test (also scaled) ###############
3
4  ## Pesaran's CD test for cross-sectional dependence in panel data models
5  ## (and Breusch and Pagan's LM and scaled LM)
6  ## ref. Pesaran, General diagnostic tests..., CESifo WP 1229, 2004
7
8  ## In case K+1>T the group-specific model is not estimable;
9  ## as in Greene 11.7.2, formula (11.23) we use the group-specific residuals
10  ## of a consistent estimator. This may be pooled OLS, RE, FE. Here the
11  ## default is set to FE.
12
13  ## Note that the test can be performed on the results of plm objects with
14  ## any kind of effects: having "time" effects means checking for
15  ## xs-dependence *after* introducing time dummies.
16
17  ## In principle, the test can be performed on the results of *any*
18  ## panelmodel object. Some issues remain regarding standardization of
19  ## model output: some missing pieces are, e.g., the 'model$indexes'
20  ## in ggls. ''fd'' models are also not compatible because of indexes
21  ## keeping the original timespan, while data lose the first period.
22
23## production version, generic and based on plm
24
25## version 11: added test = "bcsclm"
26##
27## version 10:
28## substantial optimization for speed, now fast (few seconds) on N=3000
29## all methods pass on a pseries to pcdres()
30
31## make toy example
32#dati <- data.frame(ind=rep(1:7, 4), time=rep(1:4, each=7), x=rnorm(28),
33#                   group=rep(c(1,1,2,2,2,3,3), 4))
34#pdati <- pdata.frame(dati)
35
36#' Tests of cross-section dependence for panel models
37#'
38#' Pesaran's CD or Breusch--Pagan's LM (local or global) tests for cross
39#' sectional dependence in panel models
40#'
41#' These tests are originally meant to use the residuals of separate
42#' estimation of one time--series regression for each cross-sectional
43#' unit in order to check for cross--sectional dependence (`model = NULL`).
44#' If a different model specification (`model = "within"`, `"random"`,
45#' \ldots{}) is assumed consistent, one can resort to its residuals for
46#' testing (which is common, e.g., when the time dimension's length is
47#' insufficient for estimating the heterogeneous model).
48#'
49#' If the time
50#' dimension is insufficient and `model = NULL`, the function defaults
51#' to estimation of a `within` model and issues a warning. The main
52#' argument of this function may be either a model of class
53#' `panelmodel` or a `formula` and `data frame`; in the second case,
54#' unless `model` is set to `NULL`, all usual parameters relative to
55#' the estimation of a `plm` model may be passed on. The test is
56#' compatible with any consistent `panelmodel` for the data at hand,
57#' with any specification of `effect` (except for `test = "bcsclm"` which
58#' requires a within model with either individual or two-ways effect).
59#' E.g., specifying  `effect = "time"` or `effect = "twoways"` allows
60#' to test for residual cross-sectional dependence after the introduction
61#' of time fixed effects to account for common shocks.
62#'
63#' A **local** version of either test can be computed by supplying a
64#' proximity matrix (elements coercible to `logical`) with argument
65#' `w` which provides information on whether any pair of individuals
66#' are neighbours or not. If `w` is supplied, only neighbouring pairs
67#' will be used in computing the test; else, `w` will default to
68#' `NULL` and all observations will be used. The matrix need not be
69#' binary, so commonly used "row--standardized" matrices can be
70#' employed as well. `nb` objects from \CRANpkg{spdep} must instead be
71#' transformed into matrices by \CRANpkg{spdep}'s function `nb2mat`
72#' before using.
73#'
74#' The methods implemented are suitable also for unbalanced panels.
75#'
76#' Pesaran's CD test (`test="cd"`), Breusch and Pagan's LM test
77#' (`test="lm"`), and its scaled version (`test="sclm"`) are all
78#' described in \insertCite{PESA:04;textual}{plm} (and complemented by
79#' Pesaran (2005)). The bias-corrected scaled test (`test="bcsclm"`)
80#' is due to \insertCite{BALT:FENG:KAO:12}{plm} and only valid for
81#' within models including the individual effect (it's unbalanced
82#' version uses max(Tij) for T) in the bias-correction term).
83#' \insertCite{BREU:PAGA:80;textual}{plm} is the original source for
84#' the LM test.
85#'
86#' The test on a `pseries` is the same as a test on a pooled
87#' regression model of that variable on a constant, i.e.,
88#' `pcdtest(some_pseries)` is equivalent to `pcdtest(plm(some_var ~ 1,
89#' data = some_pdata.frame, model = "pooling")` and also equivalent to
90#' `pcdtest(some_var ~ 1, data = some_data)`, where `some_var` is
91#' the variable name in the data which corresponds to `some_pseries`.
92#'
93#' @aliases pcdtest
94#' @param x an object of class `formula`, `panelmodel`, or `pseries`
95#'     (depending on the respective interface) describing the model to
96#'     be tested,
97#' @param data a `data.frame`,
98#' @param index an optional numerical index, if `NULL`, the first two
99#'     columns of the data.frame provided in argument `data` are
100#'     assumed to be the index variables; for further details see
101#'     [pdata.frame()],
102#' @param model an optional character string indicating which type of
103#'     model to estimate; if left to `NULL`, the original
104#'     heterogeneous specification of Pesaran is used,
105#' @param test the type of test statistic to be returned. One of
106#'     \itemize{ \item `"cd"` for Pesaran's CD statistic, \item `"lm"`
107#'     for Breusch and Pagan's original LM statistic, \item `"sclm"`
108#'     for the scaled version of Breusch and Pagan's LM statistic,
109#'     \item `"bcsclm"` for the bias-corrected scaled version of
110#'     Breusch and Pagan's LM statistic, \item `"rho"` for the average
111#'     correlation coefficient, \item `"absrho"` for the average
112#'     absolute correlation coefficient,}
113#' @param w either `NULL` (default) for the global tests or -- for the
114#'     local versions of the statistics -- a `n x n` `matrix`
115#'     describing proximity between individuals, with \eqn{w_ij = a}
116#'     where \eqn{a} is any number such that `as.logical(a)==TRUE`, if
117#'     \eqn{i,j} are neighbours, \eqn{0} or any number \eqn{b} such
118#'     that `as.logical(b)==FALSE` elsewhere. Only the lower
119#'     triangular part (without diagonal) of `w` after coercing by
120#'     `as.logical()` is evaluated for neighbouring information (but
121#'     `w` can be symmetric). See also **Details** and
122#'     **Examples**,
123#' @param \dots further arguments to be passed on for model estimation to `plm`,
124#'    such as `effect` or `random.method`.
125#' @return An object of class `"htest"`.
126#' @export
127#' @references
128#'
129#' \insertRef{BALT:FENG:KAO:12}{plm}
130#'
131#' \insertRef{BREU:PAGA:80}{plm}
132#'
133#' \insertRef{PESA:04}{plm}
134#'
135#' \insertRef{PESA:15}{plm}
136#'
137#' @keywords htest
138#' @examples
139#'
140#' data("Grunfeld", package = "plm")
141#' ## test on heterogeneous model (separate time series regressions)
142#' pcdtest(inv ~ value + capital, data = Grunfeld,
143#'         index = c("firm", "year"))
144#'
145#' ## test on two-way fixed effects homogeneous model
146#' pcdtest(inv ~ value + capital, data = Grunfeld, model = "within",
147#'         effect = "twoways", index = c("firm", "year"))
148#'
149#' ## test on panelmodel object
150#' g <- plm(inv ~ value + capital, data = Grunfeld, index = c("firm", "year"))
151#' pcdtest(g)
152#'
153#' ## scaled LM test
154#' pcdtest(g, test = "sclm")
155#'
156#' ## test on pseries
157#' pGrunfeld <- pdata.frame(Grunfeld)
158#' pcdtest(pGrunfeld$value)
159#'
160#' ## local test
161#' ## define neighbours for individual 2: 1, 3, 4, 5 in lower triangular matrix
162#' w <- matrix(0, ncol= 10, nrow=10)
163#' w[2,1] <- w[3,2] <- w[4,2] <- w[5,2] <- 1
164#' pcdtest(g, w = w)
165#'
166pcdtest <- function(x, ...)
167{
168    UseMethod("pcdtest")
169}
170
171## this formula method here only for adding "rho" and "absrho"
172## arguments
173
174#' @rdname pcdtest
175#' @export
176pcdtest.formula <- function(x, data, index = NULL, model = NULL,
177                            test = c("cd", "sclm", "bcsclm", "lm", "rho", "absrho"),
178                            w = NULL, ...) {
179    #data <- pdata.frame(data, index = index)
180    test <- match.arg(test)
181    if(test == "bcsclm" && (is.null(model) || model != "within"))
182      stop("for test = 'bcsclm', set argument model = 'within'")
183
184    # evaluate formula in parent frame
185    cl <- match.call(expand.dots = TRUE)
186    cl$model  <- if(test != "bcsclm") "pooling" else "within"
187      if(test == "bcsclm") {
188        # check args model and effect for test = "bcsclm"
189        if(is.null(cl$effect)) cl$effect <- "individual" # make default within model is individual within
190        eff <- isTRUE(cl$effect == "individual" || cl$effect == "twoways")
191        if(model != "within" || !eff) stop("for test = 'bcsclm', requirement is model = \"within\" and effect = \"individual\" or \"twoways\"")
192      }
193    names(cl)[2L] <- "formula"
194    m <- match(plm.arg, names(cl), 0L)
195    cl <- cl[c(1L, m)]
196    cl[[1L]] <- as.name("plm")
197    mymod <- eval(cl, parent.frame()) # mymod is either "pooling" or "within" (the latter iff for test = "bcsclm")
198
199    hetero.spec <- if(is.null(model)) TRUE else FALSE
200
201    if(hetero.spec && min(pdim(mymod)$Tint$Ti) < length(mymod$coefficients)+1) {
202      warning("Insufficient number of observations in time to estimate heterogeneous model: using within residuals",
203            call. = FALSE)
204      hetero.spec <- FALSE
205      model <- "within"
206    }
207
208    ind0 <- attr(model.frame(mymod), "index")
209    tind <- as.numeric(ind0[[2L]])
210    ind <- as.numeric(ind0[[1L]])
211
212    if(hetero.spec) {
213        ## estimate individual normal regressions one by one
214        ## (original heterogeneous specification of Pesaran)
215        X <- model.matrix(mymod)
216        y <- model.response(model.frame(mymod))
217        unind <- unique(ind)
218        n <- length(unind)
219        ti.res   <- vector("list", n)
220        ind.res  <- vector("list", n)
221        tind.res <- vector("list", n)
222        for (i in 1:n) {
223            tX <- X[ind == unind[i], , drop = FALSE]
224            ty <- y[ind == unind[i]]
225            res.i <- lm.fit(tX, ty)$residuals
226            ti.res[[i]] <- res.i
227            names(ti.res[[i]]) <- tind[ind == unind[i]]
228            ind.res[[i]] <- rep(i, length(res.i))
229            tind.res[[i]] <- tind[ind == unind[i]]
230        }
231        ## make pseries of (all) residuals
232        resdata <- data.frame(ee   = unlist(ti.res,   use.names = FALSE),
233                              ind  = unlist(ind.res,  use.names = FALSE),
234                              tind = unlist(tind.res, use.names = FALSE))
235        pee <- pdata.frame(resdata, index = c("ind", "tind"))
236        tres <- pee$ee
237    } else {
238      # else case is one of:
239      # a) insufficient number of observations for heterogen. spec. or
240      # b) model specified when function was called (incl. case test = "bcsclm")
241      if(test != "bcsclm") {
242        # Estimate the model specified originally in function call or due to
243        # forced model switch to within model by insufficient number of
244        # observations for heterogen. spec.
245        # (for test = "bcsclm" it is ensured that a within model was already
246        # estimated -> no need to estimate again a within model)
247        cl$model <- model
248        mymod <- eval(cl, parent.frame())
249      }
250
251      tres <- resid(mymod)
252      unind <- unique(ind)
253      n <- length(unind)
254      t <- min(pdim(mymod)$Tint$Ti)
255      nT <- length(ind)
256      k <- length(mymod$coefficients)
257      }
258
259    return(pcdres(tres = tres, n = n, w = w,
260                  form = paste(deparse(x)),
261                  test = test))
262}
263
264
265## panelmodel method: just fetch resid (as a pseries) and hand over to pcdres
266
267#' @rdname pcdtest
268#' @export
269pcdtest.panelmodel <- function(x, test = c("cd", "sclm", "bcsclm", "lm", "rho", "absrho"),
270                               w = NULL, ...) {
271
272    test <- match.arg(test)
273    model <- describe(x, "model")
274    effect <- describe(x, "effect")
275    eff <- (effect == "individual" || effect == "twoways")
276    if (test == "bcsclm")
277      if (model != "within" || !eff) stop("for test = 'bcsclm', model x must be a within individual or twoways model")
278
279    tres <- resid(x)
280    index <- attr(model.frame(x), "index")
281    #tind <- as.numeric(index[[2L]])
282    ind <- as.numeric(index[[1L]])
283    unind <- unique(ind)
284    n <- length(unind)
285    #t <- pdim(x)$Tint$Ti
286    #nT <- length(ind)
287    #k <- length(x$coefficients)
288    return(pcdres(tres = tres, n = n, w = w,
289                  form = paste(deparse(x$formula)),
290                  test = test))
291}
292
293#' @rdname pcdtest
294#' @export
295pcdtest.pseries <- function(x, test = c("cd", "sclm", "bcsclm", "lm", "rho", "absrho"),
296                             w = NULL, ...) {
297
298    ## calculates local or global CD test on a pseries 'x' just as it
299    ## would on model residuals
300    ## important difference here: a pseries _can_ have NAs
301
302    # input check
303    if (!inherits(x, "pseries")) stop("input 'x' needs to be of class \"pseries\"")
304    form <- paste(deparse(substitute(x)))
305
306    pos.na <- is.na(x)
307    if (any(pos.na)) {
308      x <- subset_pseries(x, !pos.na) # TODO: use [.pseries (pseries subsetting) once implemented
309      warning("NA values encountered in input and removed")
310      if (length(x) == 0L) stop("input is empty after removal of NA values")
311    }
312
313    ## get indices
314    tind <- as.numeric(attr(x, "index")[[2L]])
315    ind <- as.numeric(attr(x, "index")[[1L]])
316
317    ## det. number of groups and df
318    unind <- unique(ind)
319    n <- length(unind)
320
321    tres <- x
322
323    ## "pre-allocate" an empty list of length n
324    #tres <- vector("list", n)
325
326    ## use model residuals, group by group
327    ## list of n:
328    ## t_i residuals for each x-sect. 1..n
329    #for(i in 1:n) {
330    #          # remove NAs
331    #          xnonna <- !is.na(x[ind==unind[i]])
332    #          tres[[i]] <- x[ind==unind[i]][xnonna]
333    #          ## name resids after the time index
334    #          names(tres[[i]]) <- tind[ind==unind[i]][xnonna]
335    #          }
336
337    return(pcdres(tres = tres, n = n, w = w,
338                  form = form,
339                  test = match.arg(test)))
340}
341
342pcdres <- function(tres, n, w, form, test) {
343  # 'form' is a character describing the formula (not a formula object!)
344  # and goes into htest_object$data.name
345
346  ## Take model residuals as pseries, and calc. test
347  ## (from here on, what's needed for rho_ij is ok)
348
349  ## this function is the modulus calculating the test,
350  ## to be called from pcdtest.formula,
351  ## pcdtest.panelmodel or pcdtest.pseries
352
353  ## now (since v10) tres is the pseries of model residuals
354
355    ## calc matrix of all possible pairwise corr.
356    ## coeffs. (200x speedup from using cor())
357    wideres <- t(preshape(tres, na.rm = FALSE))
358    rho <- cor(wideres, use = "pairwise.complete.obs")
359
360    ## find length of intersecting pairs
361    ## fast method, times down 200x
362    data.res <- data.frame(time = attr(tres, "index")[[2L]],
363                           indiv = attr(tres, "index")[[1L]])
364    ## tabulate which obs in time for each ind are !na
365    presence.tab <- table(data.res)
366    ## calculate t.ij
367    t.ij <- crossprod(presence.tab)
368
369  # input check
370  if (!is.null(w)) {
371    dims.w <- dim(w)
372    if(dims.w[1L] != n || dims.w[2L] != n)
373      stop(paste0("matrix 'w' describing proximity of individuals has wrong dimensions: ",
374           "should be ", n, " x ", n, " (no. of individuals) but is ", dims.w[1L], " x ", dims.w[2L]))
375  }
376
377
378  ## begin features for local test ####################
379  ## higher orders omitted for now, use wlag() explicitly
380
381  ## if global test, set all elements in w to 1
382  if(is.null(w)) {
383    w <- matrix(1, ncol = n, nrow = n)
384    dep <- ""
385  } else { dep <- "local" }
386
387  ## make (binary) selector matrix based on the contiguity matrix w
388  ## and extracting elements corresponding to ones in the lower triangle
389  ## excluding the diagonal
390
391  ## transform in logicals (0=FALSE, else=TRUE: no need to worry
392  ## about row-std. matrices)
393  selector.mat <- matrix(as.logical(w), ncol = n)
394
395  ## some sanity checks for 'w' (not perfect sanity, but helps)
396  if (sum(selector.mat[lower.tri(selector.mat, diag = FALSE)]) == 0) {
397    stop(paste0("no neighbouring individuals defined in proximity matrix 'w'; ",
398                "only lower triangular part of 'w' (w/o diagonal) is evaluated"))
399  } else {
400    if (sum(selector.mat[upper.tri(selector.mat, diag = FALSE)]) != 0) {
401      if (!isSymmetric((unname(selector.mat)))) { # unname needed to ignore rownames and colnames
402        stop(paste0("proximity matrix 'w' is ambiguous: upper and lower triangular part ",
403                    "define different neighbours (it is sufficient to provide information ",
404                    "about neighbours only in the lower triangluar part of 'w'"))
405      }
406    }
407  }
408
409  ## if no intersection or only 1 shared period of e_it and e_jt
410  ## => exclude from calculation and issue a warning.
411  ## In general, length(m.ij) gives the number of shared periods by indiviudals i, j
412  ## Thus, non intersecting pairs are indicated by length(m.ij) == 0 (t.ij[i,j] == 0)
413  no.one.intersect <- (t.ij <= 1)
414  if (any(no.one.intersect, na.rm = TRUE)) {
415    # t.ij is a lower triangular matrix: do not divide by 2 to get the number of non-intersecting pairs!
416    number.of.non.one.intersecting.pairs <- sum(no.one.intersect, na.rm = TRUE)
417    number.of.total.pairs <- (n*(n-1))/2
418    share.on.one.intersect.pairs <- number.of.non.one.intersecting.pairs / number.of.total.pairs * 100
419    warning(paste("Some pairs of individuals (",
420                  signif(share.on.one.intersect.pairs, digits = 2),
421                  " percent) do not have any or just one time period in common and have been omitted from calculation", sep=""))
422    selector.mat[no.one.intersect] <- FALSE
423  }
424
425  ## set upper tri and diagonal to FALSE
426  selector.mat[upper.tri(selector.mat, diag = TRUE)] <- FALSE
427
428  ## number of elements in selector.mat
429  ## elem.num = 2*(N*(N-1)) in Pesaran (2004), formulae (6), (7), (31), ...
430  elem.num <- sum(selector.mat)
431
432  ## end features for local test ######################
433
434  ## Breusch-Pagan or Pesaran statistic for cross-sectional dependence,
435  ## robust vs. unbalanced panels:
436
437  switch(test,
438   lm = {
439    CDstat        <- sum((t.ij*rho^2)[selector.mat])
440    pCD           <- pchisq(CDstat, df = elem.num, lower.tail = FALSE)
441    names(CDstat) <- "chisq"
442    parm          <- elem.num
443    names(parm)   <- "df"
444    testname      <- "Breusch-Pagan LM test"
445   },
446   sclm = {
447    CDstat        <- sqrt(1/(2*elem.num))*sum((t.ij*rho^2-1)[selector.mat])
448    pCD           <- 2*pnorm(abs(CDstat), lower.tail = FALSE)
449    names(CDstat) <- "z"
450    parm          <- NULL
451    testname      <- "Scaled LM test"
452   },
453   bcsclm = {
454     # Baltagi/Feng/Kao (2012), formula (11)
455     # (unbalanced case as sclm + in bias correction as EViews: max(T_ij) instead of T)
456      CDstat        <- sqrt(1/(2*elem.num))*sum((t.ij*rho^2-1)[selector.mat]) - (n/(2*(max(t.ij)-1)))
457      pCD           <- 2*pnorm(abs(CDstat), lower.tail = FALSE)
458      names(CDstat) <- "z"
459      parm          <- NULL
460      testname      <- "Bias-corrected Scaled LM test"
461   },
462   cd = {
463     # (Pesaran (2004), formula (31))
464    CDstat        <- sqrt(1/elem.num)*sum((sqrt(t.ij)*rho)[selector.mat])
465    pCD           <- 2*pnorm(abs(CDstat), lower.tail = FALSE)
466    names(CDstat) <- "z"
467    parm          <- NULL
468    testname      <- "Pesaran CD test"
469   },
470   rho = {
471    CDstat        <- sum(rho[selector.mat])/elem.num
472    pCD           <- NULL
473    names(CDstat) <- "rho"
474    parm          <- NULL
475    testname      <- "Average correlation coefficient"
476   },
477   absrho = {
478    CDstat        <- sum(abs(rho)[selector.mat])/elem.num
479    pCD           <- NULL
480    names(CDstat) <- "|rho|"
481    parm          <- NULL
482    testname      <- "Average absolute correlation coefficient"
483   })
484
485  ##(insert usual htest features)
486  RVAL <- list(statistic = CDstat,
487               parameter = parm,
488               method    = paste(testname, "for", dep,
489                            "cross-sectional dependence in panels"),
490               alternative = "cross-sectional dependence",
491               p.value     = pCD,
492               data.name   = form)
493  class(RVAL) <- "htest"
494  return(RVAL)
495}
496
497preshape <- function(x, na.rm = TRUE, ...) {
498    ## reshapes pseries,
499    ## e.g., of residuals from a panelmodel,
500    ## in wide form
501    inames <- names(attr(x, "index"))
502    mres <- reshape(cbind(as.vector(x), attr(x, "index")),
503                    direction = "wide",
504                    timevar = inames[2L],
505                    idvar = inames[1L])
506    ## drop ind in first column
507    mres <- mres[ , -1L, drop = FALSE]
508    ## reorder columns (may be scrambled depending on first
509    ## available obs in unbalanced panels)
510    mres <- mres[ , order(dimnames(mres)[[2L]])]
511    ## if requested, drop columns (time periods) with NAs
512    if(na.rm) {
513        na.cols <- vapply(mres, FUN = anyNA, FUN.VALUE = TRUE, USE.NAMES = FALSE)
514        if(sum(na.cols) > 0L) mres <- mres[ , !na.cols]
515    }
516    return(mres)
517}
518
519
520
521
522#' Cross--sectional correlation matrix
523#'
524#' Computes the cross--sectional correlation matrix
525#'
526#'
527#' @param x an object of class `pseries`
528#' @param grouping grouping variable,
529#' @param groupnames a character vector of group names,
530#' @param value to complete,
531#' @param \dots further arguments.
532#' @return A matrix with average correlation coefficients within a group
533#' (diagonal) and between groups (off-diagonal).
534#' @export
535#' @keywords htest
536#' @examples
537#'
538#' data("Grunfeld", package = "plm")
539#' pGrunfeld <- pdata.frame(Grunfeld)
540#' grp <- c(rep(1, 100), rep(2, 50), rep(3, 50)) # make 3 groups
541#' cortab(pGrunfeld$value, grouping = grp, groupnames = c("A", "B", "C"))
542#'
543cortab <- function(x, grouping, groupnames = NULL,
544                   value = "statistic", ...) {
545    ## makes matrix containing within (diagonal) and between (off-diagonal)
546    ## correlation
547    ## needs a pseries and a groupings vector of **same length**
548
549    ## would use a better naming, and also passing a char or factor as
550    ## grouping index
551
552    ## x must be a pseries
553    if(!inherits(x, "pseries")) stop("First argument must be a pseries")
554    if(length(x) != length(grouping)) stop("arguments 'x' and 'grouping' must have same length")
555
556    fullind <- as.numeric(attr(x, "index")[ , 1L])
557    ids <- unique(fullind)
558    n <- length(ids)
559    regs <- 1:length(unique(grouping))
560
561    if(!(is.numeric(grouping))) grouping <- as.numeric(as.factor(grouping))
562
563    idnames <- as.character(ids)
564    if(is.null(groupnames)) {
565        groupnames <- as.character(unique(grouping))
566    }
567
568    ## make matrices of between-regions correlations
569    ## (includes within correlation on diagonal)
570    ## for each pair of regions (nb: no duplicates, e.g., 3.1 but not 1.3)
571
572    ## make w<1.n>:
573    for(h in 1:length(regs)) {
574      for(k in 1:h) {
575        statew <- matrix(0, ncol = n, nrow = n)
576        ## make statew for cor. between h and k
577        for(i in 1:n) {
578          ## get first region (all values equal, so take first one)
579          ireg <- grouping[fullind == ids[i]][1L]
580          if(ireg == h) {
581            for(j in 1:n) {
582                jreg <- grouping[fullind == ids[j]][1L]
583                if(jreg == k) statew[i, j] <- 1
584            }
585          }
586        }
587        if(h!=k) statew <- statew + t(statew)
588        ## just for debugging reasons:
589        dimnames(statew) <- list(idnames, idnames)
590        ## eliminate self.correlation of states if i=j
591        diag(statew) <- 0
592        ## not needed: pcdtest seems to do this by construction
593        eval(parse(text=paste("w", h, ".", k, " <- statew", sep="")))
594      }
595     }
596
597     ## notice: without the line
598     ## '' if(i!=j) statew <- statew + t(statew) ''
599     ## all wn.n matrices would have values only on one half (upper
600     ## or lower triangle)
601
602     ## make generic table of regions' within and between correlation
603     ## argument: a pseries
604    #YC regnames is undefined, so is myw
605    tab.g <- function(x, regs, regnames, test="rho", value) {
606        myw <- 0
607         tabg <- matrix(NA, ncol=length(regs), nrow=length(regs))
608         for(i in 1:length(regs)) {
609             for(j in 1:i) {
610                 ## take appropriate w matrix
611                 eval(parse(text = paste("myw<-w", i, ".", j, sep = "")))
612                 tabg[i, j] <- pcdtest(x, test = "rho", w = myw)[[value]]
613             }
614         }
615         dimnames(tabg) <- list(groupnames, groupnames)
616         return(tabg)
617    }
618    regnames <- ""
619    mytab <- tab.g(x, regs = regs, regnames = regnames, test = "rho", value = value)
620    return(mytab)
621}
622
623
624
625
626