diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 46c7d20..77eaa7e 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -15,6 +15,10 @@ jobs: env: R_KEEP_PKG_SOURCE: yes R_PROFILE_USER: /dev/null + # pak workers read the project profile; CI uses the action's library. + RENV_CONFIG_AUTOLOADER_ENABLED: "false" + # Allow slower Windows workers up to 60 seconds to start (milliseconds). + PKG_SUBPROCESS_TIMEOUT: "60000" steps: - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/NAMESPACE b/NAMESPACE index c48de94..ffe4fe1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,7 +2,9 @@ S3method(augment,csdm_fit) S3method(cd_test,csdm_fit) +S3method(cd_test,data.frame) S3method(cd_test,default) +S3method(cd_test,pdata.frame) S3method(coef,csdm_fit) S3method(df.residual,csdm_fit) S3method(fitted,csdm_fit) @@ -13,6 +15,7 @@ S3method(model.matrix,csdm_fit) S3method(nobs,csdm_fit) S3method(predict,csdm_fit) S3method(print,cd_test) +S3method(print,cd_test_list) S3method(print,csdm_fit) S3method(print,summary.csdm_fit) S3method(residuals,csdm_fit) diff --git a/NEWS.md b/NEWS.md index 03572ea..9dd53be 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,24 @@ +# csdm (development version) + +- Add `cd_test()` methods for `data.frame` and `pdata.frame` inputs. Select + numeric panel variables explicitly using bare names, quoted names, or + character vectors. Plain data frames require named `id` and `time` columns; + indexed panels use their stored indexes. Return separate named results per + variable and print a combined table with variable-specific sample dimensions. +- Add `reps` to CDw and CDw+, defaulting to 1. Repeated weighted statistics use + the sum divided by `sqrt(reps)`, following equation (33) of Juodis and Reese; + CDw+ adds its screening term once. Seeded calls continue to preserve RNG state. +- Record incomplete time periods removed by `drop.incomplete.times` in + `excluded_times`, alongside periods with no finite observations. +- Correct fitted-model CD* input to retain the common-factor component after + subtracting the economic and deterministic fitted terms, following + Pesaran-Xie's regression input construction. This changes fitted CD* results; + CD, CDw, and CDw+ still use full fitted residuals. Preserve and document the + standardized-PCA implementation variant for CD*. +- Clarify the CDw and CDw+ formulas against Juodis-Reese equations (30) and (32), + document the complete CD* correction, and correct the Pesaran-Xie and + Fan-Liao-Yao references. + # csdm 2.0.0 Version 2.0.0 expands the package's R model interfaces, adds new sample diff --git a/R/utils_cd.R b/R/utils_cd.R index 3bd6b7f..cab7c04 100644 --- a/R/utils_cd.R +++ b/R/utils_cd.R @@ -1,53 +1,115 @@ # utils_cd.R - Cross-sectional dependence tests for panel models -#' Cross-sectional dependence (CD) tests for panel residuals +#' Cross-sectional dependence (CD) tests for panel data and residuals #' #' Computes Pesaran CD, CDw, CDw+, and CD* tests for cross-sectional dependence -#' in panel residuals. The implementation supports residual matrices or fitted -#' \code{csdm_fit} objects and provides consistent handling of unbalanced panels. +#' in panel variables or residuals. The implementation supports data frames, +#' indexed panel data frames, numeric matrices, and fitted \code{csdm_fit} objects. #' -#' @param object A \code{csdm_fit} model object or a numeric matrix of residuals (N x T). -#' @param ... Additional arguments passed to methods. +#' @param object A \code{data.frame}, \code{pdata.frame}, \code{csdm_fit} model, +#' or numeric matrix with units in rows and time periods in columns (N x T). +#' @param ... For data frames, explicitly selected numeric columns as bare names, +#' quoted names, or character vectors of names. Selections must be unnamed. +#' For matrices and fitted models, additional arguments passed to methods. #' #' @return An object of class \code{cd_test} with fields \code{tests}, \code{type}, #' \code{N}, \code{T}, \code{na.action}, \code{excluded_units}, #' \code{excluded_times}, \code{kept_times}, and \code{call}. The \code{tests} #' list contains one or more test results, each with \code{statistic} and #' \code{p.value}. +#' Data-frame methods return a \code{cd_test_list}: a named list containing one +#' \code{cd_test} result per selected variable, including when only one variable +#' is selected. Each result also records \code{variable}, \code{units}, and +#' \code{excluded_unit_ids}. The list has \code{call}, \code{id}, and \code{time} +#' attributes and prints a combined table with each variable's sample dimensions. +#' For fitted models, CD* additionally records \code{input = "partial_residuals"} +#' and its own sample exclusions in \code{tests$CDstar}. With \code{type = "all"}, +#' top-level sample fields describe the full residual panel; the CD* sample +#' dimensions are recorded in \code{tests$CDstar$N} and \code{tests$CDstar$T}. #' #' @details +#' ## Selecting panel variables +#' +#' For a plain data frame, supply distinct unit and time column names through +#' \code{id} and \code{time}. For a \code{pdata.frame}, the stored indexes are +#' used, even when the index columns have been removed from the data. Select at +#' least one numeric non-index column explicitly. Each variable is tested +#' separately using its own available sample and the same missing-data policy. +#' No regression is fitted. CDw demeans the observations within each unit, as it +#' does for matrix inputs. Test controls and panel indexes follow \code{...} in +#' the data-frame methods and must be named. +#' #' ## Notation #' -#' Let \eqn{E} be the residual matrix with \eqn{N} cross-sectional units and \eqn{T} -#' time periods. For each unit pair \eqn{(i,j)}, let \eqn{T_{ij}} be the number of -#' overlapping time periods and \eqn{\rho_{ij}} the pairwise correlation. +#' Let \eqn{E} contain the selected observations or residuals, with \eqn{N} +#' cross-sectional units and \eqn{T} time periods. For each unit pair +#' \eqn{(i,j)}, let \eqn{T_{ij}} be the number of +#' overlapping time periods and \eqn{\widehat\rho_{ij}} the correlation computed +#' after demeaning both series over that pair's overlapping observations. #' #' ## Test statistics #' #' \describe{ #' \item{CD (Pesaran, 2015)}{ -#' \deqn{CD = \sqrt{\frac{2}{N(N-1)}} \sum_{i 2\sqrt{\log(N)/T}\right\}.} -#' Here \eqn{\rho_{ij}} is the ordinary residual correlation, without +#' Equation (32) of Juodis and Reese applies the power-enhancement principle +#' of Fan, Liao, and Yao. Its statistic is +#' \deqn{CD_{W+} = CD_W + \sum_{i 2\sqrt{\log(N)/T}\right\}.} +#' Here \eqn{\widehat\rho_{ij}} is the ordinary residual correlation, without #' multiplication by \eqn{\sqrt{T}}. The nonnegative screening term is #' asymptotically zero under the conditions of the null hypothesis. +#' With \code{reps > 1}, replace \eqn{CD_W} by \eqn{\overline{CD}_W} and add +#' the screening term once. #' } -#' \item{CD* (Pesaran and Xie, 2021)}{ -#' CD is computed on residuals after removing \code{n_pc} principal components -#' from \eqn{E}. This provides a bias-corrected test under multifactor errors. +#' \item{CD* (Pesaran-Xie correction; standardized-PCA variant)}{ +#' Each unit is demeaned and divided by its sample standard deviation before +#' extracting \code{n_pc} principal components. Write \eqn{A} for this +#' \eqn{T\times N} standardized matrix. The estimated loading matrix +#' \eqn{\widehat\Gamma} is normalized so that +#' \eqn{\widehat\Gamma'\widehat\Gamma/N=I}, and the factors are +#' \eqn{\widehat F=A\widehat\Gamma/N}. For the factor-filtered residuals +#' \eqn{r_{it}} of \eqn{A-\widehat F\widehat\Gamma'}, define +#' \deqn{\widehat\sigma_i=\sqrt{T^{-1}\sum_t r_{it}^2},\quad +#' \widehat\varphi=N^{-1}\sum_i\widehat\gamma_i/\widehat\sigma_i,\quad +#' \widehat a_i=1-\widehat\sigma_i\widehat\varphi'\widehat\gamma_i.} +#' The implemented bias correction is +#' \deqn{\widehat\theta=1-N^{-1}\sum_i\widehat a_i^2,\qquad +#' CD^*=\frac{CD(r)+\sqrt{T/2}\,\widehat\theta}{1-\widehat\theta}.} +#' This uses the Pesaran-Xie correction algebra, but their PCA procedure +#' operates on unstandardized observations. Standardizing before PCA can +#' change the estimated factor space and is an implementation variant. +#' Setting \code{n_pc = 0} returns classical CD without factor removal. +#' For fitted models, the PCA input is +#' \eqn{\widehat v_{it}=y_{it}-x_{it}'\widehat\beta_i}, where \eqn{x_{it}} +#' includes all economic and deterministic design columns, including any +#' intercept, trend, and constructed lags. The fitted CSA contribution is +#' retained, following the input construction in Pesaran-Xie's regression +#' procedure. CD, CDw, and CDw+ use the full regression residuals instead. #' } #' } #' @@ -56,6 +118,16 @@ #' loading/error-scale designs after standardization. Numerical rank checks #' do not establish the validity of the asymptotic approximation. #' +#' All p-values use the two-sided standard-normal approximation +#' \eqn{2\Phi(-|\mathrm{statistic}|)}. The implemented tests do not estimate a +#' serial-correlation correction. Stationarity, time-series dependence, factor +#' strength, and relative panel dimensions must satisfy the relevant theory; +#' accepting a numeric variable does not establish these conditions. For CD, +#' omitting pairs or using very short overlaps can also affect calibration. +#' Classical CD can be biased on CCE residuals because common time parameters +#' have been estimated. CDw addresses this problem under Juodis-Reese's +#' assumptions, including \eqn{\sqrt{T}/N\to0}. +#' #' ## Missing data and balance #' #' Time periods containing no finite residual for any retained unit are outside @@ -82,6 +154,8 @@ #' #' \insertRef{PesaranXie2021}{csdm} #' +#' \insertRef{PesaranXie2026}{csdm} +#' #' @examples #' # Simulate independent and dependent panels #' set.seed(1) @@ -95,6 +169,12 @@ #' # Specific test with parameters #' cd_test(E_indep, type = "CDstar", n_pc = 2) #' +#' # Test raw panel variables separately +#' panel <- expand.grid(id = 1:10, year = 1:10) +#' panel$x <- rnorm(nrow(panel)) +#' panel$y <- rnorm(nrow(panel)) +#' cd_test(panel, x, "y", id = "id", time = "year") +#' #' # From a fitted csdm model #' data(PWT_60_07, package = "csdm") #' df <- PWT_60_07 @@ -120,8 +200,14 @@ cd_test <- function(object, ...) { #' \code{"CDstar"}, or \code{"all"} (default: \code{"CD"}). #' @param n_pc Number of principal components for CD* (default 4). #' @param seed Integer seed for weight draws. Seeded calls restore the caller's RNG state; NULL uses the current RNG stream. -#' @param min_overlap Minimum number of overlapping time periods required for a unit -#' pair to be included in CD/CDw/CDw+ (default 2). +#' @param reps Positive integer number of independent weight draws for CDw and +#' CDw+ (default 1). The draws are summed and divided by \code{sqrt(reps)}. +#' CD and CD* are unaffected. Weighted test results record the number of draws +#' in their \code{reps} field. For data inputs, \code{seed} is applied separately +#' to each variable's test, making seeded results independent of selection order. +#' @param min_overlap Minimum number of finite observations for retaining a unit +#' before time filtering (all tests), and minimum overlapping time periods +#' for including a pair in classical CD (default 2). #' @param na.action How to handle missing data: \code{"drop.incomplete.times"} #' removes time periods with any missing observations to create a balanced panel for CD*; #' \code{"pairwise"} (default) uses pairwise correlations for CD; unbalanced CDw/CDw+ requests error and CD* warns. @@ -133,7 +219,22 @@ cd_test.default <- function(object, seed = NULL, min_overlap = 2L, na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L, ...) { + result <- .cd_test_matrix(object, type = type, n_pc = n_pc, seed = seed, + min_overlap = min_overlap, na.action = na.action, + reps = reps, ...) + result$call <- match.call() + result +} + +# The fitted-model method supplies a different panel to CD*, so it can request +# the other tests without first removing factors from the full CCE residuals. +.cd_test_matrix <- function(object, + type = c("CD", "CDw", "CDw+", "CDstar", "all"), + n_pc = 4L, seed = NULL, min_overlap = 2L, + na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L, .compute_star = TRUE, ...) { type <- match.arg(type) na.action <- match.arg(na.action) @@ -147,7 +248,9 @@ cd_test.default <- function(object, min_overlap <- .csdm_integer(min_overlap, "min_overlap") if (min_overlap < 2L) stop("'min_overlap' must be at least two.") - if (any(is.infinite(object))) stop("Residuals may contain NA, but not infinite values.") + reps <- .csdm_integer(reps, "reps") + if (reps < 1L) stop("'reps' must be at least one.", call. = FALSE) + if (any(is.infinite(object))) stop("Input values may contain NA, but not infinite values.") if (!is.null(seed)) { seed <- .csdm_integer(seed, "seed") had_seed <- exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE) @@ -165,9 +268,11 @@ cd_test.default <- function(object, excluded_units <- which(!usable) E <- object[usable, , drop = FALSE] - # Periods with no estimated residuals are not part of the residual sample. + # Periods with no finite values are not part of the diagnostic sample. + time_labels <- if (is.null(colnames(E))) seq_len(ncol(E)) else colnames(E) empty_times <- colSums(is.finite(E)) == 0L - excluded_times <- if (is.null(colnames(E))) which(empty_times) else colnames(E)[empty_times] + excluded_times <- time_labels[empty_times] + kept_time_indices <- which(!empty_times) if (any(empty_times)) E <- E[, !empty_times, drop = FALSE] if (ncol(E) < 2L) stop("At least two time periods with residual observations are required.") @@ -179,6 +284,9 @@ cd_test.default <- function(object, if (sum(complete_times) < 2) { stop("cd_test: After dropping incomplete time periods, fewer than 2 periods remain.") } + excluded <- empty_times + excluded[kept_time_indices[!complete_times]] <- TRUE + excluded_times <- time_labels[excluded] E <- E[, complete_times, drop = FALSE] if (n_dropped > 0) { message(sprintf("cd_test: Dropped %d incomplete time period%s (%.1f%%). Balanced panel: %d units x %d periods.", @@ -210,15 +318,18 @@ cd_test.default <- function(object, if (type %in% c("CDw", "CDw+", "all")) { if (anyNA(E)) stop("CDw/CDw+ currently require a balanced sample; use explicit drop.incomplete.times or classical CD.") if (!is.null(seed)) set.seed(seed) - w <- sample(c(-1, 1), N, replace = TRUE) centered <- sweep(data_tn, 2L, colMeans(data_tn)) variance <- mean(centered^2) - weighted <- sweep(centered, 2L, w, "*") - pair_sum <- sum(rowSums(weighted)^2 - rowSums(weighted^2)) / 2 - statistic <- sqrt(2 / (Tt * N * (N - 1))) * pair_sum / variance + draws <- vapply(seq_len(reps), function(g) { + w <- sample(c(-1, 1), N, replace = TRUE) + weighted <- sweep(centered, 2L, w, "*") + pair_sum <- sum(rowSums(weighted)^2 - rowSums(weighted^2)) / 2 + sqrt(2 / (Tt * N * (N - 1))) * pair_sum / variance + }, numeric(1)) + statistic <- sum(draws) / sqrt(reps) out$CDw <- list(statistic = statistic, p.value = 2 * stats::pnorm(abs(statistic), lower.tail = FALSE), - N = N, T = Tt, pairs_used = choose(N, 2)) + N = N, T = Tt, pairs_used = choose(N, 2), reps = reps) if (type %in% c("CDw+", "all")) { correlation <- stats::cor(centered) rho <- abs(correlation[upper.tri(correlation)]) @@ -227,13 +338,13 @@ cd_test.default <- function(object, statistic <- statistic + enhancement out$CDw_plus <- list(statistic = statistic, p.value = 2 * stats::pnorm(abs(statistic), lower.tail = FALSE), - N = N, T = Tt, pairs_used = choose(N, 2), + N = N, T = Tt, pairs_used = choose(N, 2), reps = reps, threshold = threshold, enhancement = enhancement) } } # 4. CD* (bias-corrected with PCA factor removal) - if (type %in% c("CDstar", "all")) { + if (.compute_star && type %in% c("CDstar", "all")) { # CD* requires balanced panel if (anyNA(E)) { warning("cd_test: CD* requires a balanced panel; returning NA due to missing values.") @@ -280,16 +391,70 @@ cd_test.csdm_fit <- function(object, seed = NULL, min_overlap = 2L, na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L, ...) { type <- match.arg(type) na.action <- match.arg(na.action) E <- .csdm_get_residuals(object, type = "auto", strict = TRUE) - cd_test.default(E, type = type, n_pc = n_pc, seed = seed, - min_overlap = min_overlap, na.action = na.action, ...) + if (type %in% c("CDstar", "all")) { + V <- .cdstar_fit_input(object, E) + star <- cd_test.default(V, type = "CDstar", n_pc = n_pc, seed = seed, + min_overlap = min_overlap, na.action = na.action, + reps = reps, ...) + star$tests$CDstar$input <- "partial_residuals" + # CD* can retain different units from tests on the full residuals. + star$tests$CDstar$excluded_units <- star$excluded_units + star$tests$CDstar$excluded_times <- star$excluded_times + star$tests$CDstar$kept_times <- star$kept_times + if (type == "CDstar") { + star$call <- match.call() + return(star) + } + } + result <- .cd_test_matrix(E, type = type, n_pc = n_pc, seed = seed, + min_overlap = min_overlap, na.action = na.action, + reps = reps, .compute_star = FALSE, ...) + if (type == "all") { + result$tests$CDstar <- star$tests$CDstar + result$type <- names(result$tests) + } + result$call <- match.call() + result +} + +# Pesaran-Xie start PCA from y minus the economic/deterministic component, +# retaining the common-factor component rather than using full CCE residuals. +.cdstar_fit_input <- function(object, E) { + X <- object$model_matrix + y <- if (!is.null(object$model_frame)) stats::model.response(object$model_frame) + rows <- object$sample$row[object$sample$used] + if (!is.matrix(X) || !is.numeric(y) || nrow(X) != length(rows) || + length(y) != length(rows) || !is.matrix(object$coef_i)) { + stop("CD* requires stored economic design and unit coefficients; refit the model.") + } + keys <- object$data[rows, c(object$id, object$time), drop = FALSE] + ids <- as.character(keys[[object$id]]) + partial <- numeric(length(rows)) + for (uid in unique(ids)) { + selected <- which(ids == uid) + beta <- object$coef_i[uid, colnames(X)] + partial[selected] <- y[selected] - as.numeric(X[selected, , drop = FALSE] %*% beta) + } + V <- E + V[] <- NA_real_ + times <- keys[[object$time]] + if (isTRUE(object$meta$pdata) && is.factor(times)) { + times <- suppressWarnings(as.numeric(as.character(times))) + } + cells <- cbind(match(ids, rownames(V)), + match(as.character(times), colnames(V))) + if (anyNA(cells)) stop("Stored CD* sample indexes are inconsistent; refit the model.") + V[cells] <- partial + V } #' @rdname cd_test -#' @param x An object of class \code{cd_test}. +#' @param x An object of class \code{cd_test} or \code{cd_test_list}. #' @param digits Number of digits to print (default 3). #' @export #' @method print cd_test @@ -300,7 +465,10 @@ print.cd_test <- function(x, digits = 3, ...) { } cat("Cross-sectional dependence tests\n") - if (!is.null(x$N) && !is.null(x$T)) { + different_samples <- any(vapply(x$tests, function(test) { + !is.null(test$N) && !is.null(test$T) && (test$N != x$N || test$T != x$T) + }, logical(1))) + if (!different_samples && !is.null(x$N) && !is.null(x$T)) { cat(sprintf("N = %d, T = %d\n", x$N, x$T)) } cat("\n") @@ -332,6 +500,10 @@ print.cd_test <- function(x, digits = 3, ...) { row.names = display_names, stringsAsFactors = FALSE ) + if (different_samples) { + out$N <- vapply(tests, get_num, key = "N", FUN.VALUE = numeric(1)) + out$T <- vapply(tests, get_num, key = "T", FUN.VALUE = numeric(1)) + } } else { key <- type_map[[as.character(x$type[1])]] if (is.null(key) || is.null(tests[[key]])) { diff --git a/R/utils_cd_data.R b/R/utils_cd_data.R new file mode 100644 index 0000000..111eed8 --- /dev/null +++ b/R/utils_cd_data.R @@ -0,0 +1,132 @@ +#' @rdname cd_test +#' @param id,time Column names (strings) for the unit and time indexes of a plain +#' data frame. For a \code{pdata.frame}, these are inferred from its stored +#' indexes and supplied values are ignored. Time labels need not be numeric. +#' @export +#' @method cd_test data.frame +cd_test.data.frame <- function(object, ..., id = NULL, time = NULL, + type = c("CD", "CDw", "CDw+", "CDstar", "all"), + n_pc = 4L, seed = NULL, min_overlap = 2L, + na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L) { + vars <- .cd_select_variables(as.list(substitute(list(...)))[-1L], + names(object), parent.frame()) + panel <- .cd_panel_input(object, id, time, vars) + type <- match.arg(type) + na.action <- match.arg(na.action) + call <- match.call() + results <- lapply(vars, function(variable) { + E <- matrix(NA_real_, length(panel$units), length(panel$times), + dimnames = list(panel$units, as.character(panel$times))) + E[panel$cells] <- as.numeric(object[[variable]]) + result <- cd_test.default(E, type = type, n_pc = n_pc, seed = seed, + min_overlap = min_overlap, na.action = na.action, reps = reps) + result$variable <- variable + result$units <- panel$units[!seq_along(panel$units) %in% result$excluded_units] + result$excluded_unit_ids <- panel$units[result$excluded_units] + result$call <- call + result + }) + names(results) <- vars + structure(results, class = c("cd_test_list", "list"), call = call, + id = panel$id, time = panel$time) +} + +# The data-frame method reads the pdata.frame indexes without requiring plm. +#' @rdname cd_test +#' @export +#' @method cd_test pdata.frame +cd_test.pdata.frame <- cd_test.data.frame + +.cd_select_variables <- function(expressions, columns, envir) { + if (!length(expressions)) { + stop("Select at least one variable explicitly in '...'.", call. = FALSE) + } + if (!is.null(names(expressions)) && any(nzchar(names(expressions)))) { + stop("Variable selections in '...' must be unnamed; name test controls explicitly.", + call. = FALSE) + } + selections <- lapply(expressions, function(expr) { + if (is.symbol(expr) && as.character(expr) %in% columns) { + return(as.character(expr)) + } + value <- tryCatch(eval(expr, envir), error = function(e) NULL) + if (!is.character(value) || !length(value) || anyNA(value) || any(!nzchar(value))) { + stop("Select variables using bare column names, quoted names, or character vectors.", + call. = FALSE) + } + value + }) + vars <- unlist(selections, use.names = FALSE) + if (anyDuplicated(vars)) stop("Select each variable only once.", call. = FALSE) + unknown <- setdiff(vars, columns) + if (length(unknown)) stop("Unknown variable(s): ", paste(unknown, collapse = ", "), + call. = FALSE) + vars +} + +.cd_panel_input <- function(object, id, time, vars) { + if (!nrow(object) || anyDuplicated(names(object))) { + stop("Supply nonempty data with unique column names.", call. = FALSE) + } + if (inherits(object, "pdata.frame")) { + index <- attr(object, "index") + if (!is.data.frame(index) || ncol(index) < 2L || nrow(index) != nrow(object)) { + stop("The pdata.frame must have valid stored unit and time indexes.", call. = FALSE) + } + id <- names(index)[1L] + time <- names(index)[2L] + keys <- index[1:2] + } else { + if (!is.character(id) || length(id) != 1L || is.na(id) || !nzchar(id) || + !is.character(time) || length(time) != 1L || is.na(time) || !nzchar(time) || + id == time || !all(c(id, time) %in% names(object))) { + stop("'id' and 'time' must name distinct columns in the data frame.", call. = FALSE) + } + keys <- object[c(id, time)] + } + valid_key <- function(x) { + is.atomic(x) && is.null(dim(x)) && !is.complex(x) && + !anyNA(x) && all(nzchar(as.character(x))) && + (!(is.numeric(x) || inherits(x, "Date") || inherits(x, "POSIXt")) || + all(is.finite(x))) + } + if (!all(vapply(keys, valid_key, logical(1)))) { + stop("Panel keys must be nonmissing and finite, with no empty labels.", call. = FALSE) + } + if (anyDuplicated(keys)) stop("Duplicate id/time cells are not allowed.", call. = FALSE) + if (any(vars %in% c(id, time))) { + stop("Select test variables other than the unit and time indexes.", call. = FALSE) + } + numeric_column <- function(x) is.numeric(x) && is.null(dim(x)) + if (!all(vapply(object[vars], numeric_column, logical(1)))) { + stop("Selected test variables must be numeric vectors.", call. = FALSE) + } + units <- sort(unique(as.character(keys[[1L]]))) + times <- sort(unique(keys[[2L]])) + list(id = id, time = time, units = units, times = times, + cells = cbind(match(as.character(keys[[1L]]), units), match(keys[[2L]], times))) +} + +#' @rdname cd_test +#' @export +#' @method print cd_test_list +print.cd_test_list <- function(x, digits = 3, ...) { + rows <- lapply(names(x), function(variable) { + result <- x[[variable]] + selected <- if (identical(result$type, "CDw+")) "CDw_plus" else result$type + tests <- result$tests[selected] + data.frame(variable = variable, test = gsub("CDw_plus", "CDw+", names(tests), fixed = TRUE), + statistic = vapply(tests, function(test) test$statistic, numeric(1)), + p.value = vapply(tests, function(test) test$p.value, numeric(1)), + N = result$N, T = result$T, row.names = NULL) + }) + cat("Cross-sectional dependence tests by variable\n\n") + table <- do.call(rbind, rows) + for (column in c("statistic", "p.value")) { + table[[column]] <- ifelse(is.na(table[[column]]), "NA", + formatC(table[[column]], digits = digits, format = "f")) + } + print(table, row.names = FALSE, right = TRUE) + invisible(x) +} diff --git a/R/utils_residuals.R b/R/utils_residuals.R index 148432c..b1e01d2 100644 --- a/R/utils_residuals.R +++ b/R/utils_residuals.R @@ -172,8 +172,9 @@ get_residuals <- function(object, #' @details #' This helper is retained temporarily for compatibility. Row scaling changes #' the weighted covariance underlying CDw, while cross-sectional time demeaning -#' can mechanically induce dependence. Neither transformation is applied by -#' `cd_test()`. +#' can mechanically induce dependence. CDw does not standardize unit variances +#' or demean across units. CD* uses within-unit standardization before PCA, +#' as described in [cd_test()]. #' #' ## Transformation steps #' @@ -185,9 +186,9 @@ get_residuals <- function(object, #' #' ## Why this preprocessing matters #' -#' CD-type tests are sensitive to scale heterogeneity and sparse columns in -#' unbalanced panels. This helper creates a better-conditioned input matrix while -#' preserving as much usable information as possible. +#' These transformations change the diagnostic input. Their suitability depends +#' on the selected test and its assumptions; preprocessing does not establish +#' the validity of the test's reference distribution. #' #' @examples #' data(PWT_60_07, package = "csdm") diff --git a/README.md b/README.md index 4aa9821..4f8a993 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ `csdm` estimates heterogeneous panel models when units may share unobserved common factors. It provides mean-group (MG), common correlated effects (CCE), dynamic CCE (DCCE), and cross-sectionally augmented ARDL (CS-ARDL) estimators, -along with residual cross-sectional dependence diagnostics. +along with cross-sectional dependence diagnostics for residuals and panel variables. The package follows the econometric structure used by Stata's `xtdcce2`, while using standard R model methods and explicit specification objects. It does not @@ -92,13 +92,27 @@ preserves gaps in calendar time. ## Cross-sectional dependence diagnostics -`cd_test()` accepts a fitted `csdm` model or an `N` by `T` residual matrix: +`cd_test()` accepts a fitted `csdm` model, an `N` by `T` numeric matrix, +or explicitly selected variables in a `data.frame` or `pdata.frame`: ```r cd_test(cce, type = "CD") cd_test(cce, type = "all", seed = 42) + +cd_test(dat, log_rgdpo, "log_hc", id = "id", time = "year") + +vars <- c("log_rgdpo", "log_hc") +cd_test(dat, vars, id = "id", time = "year", type = "CDw", + reps = 30, seed = 42) ``` +Data inputs return a named list of separate results per variable, printed in a +combined table. Variables must be selected explicitly through `...`, using bare +names, quoted names, or character vectors. Name `id`, `time`, and all test +controls; a `pdata.frame` supplies its stored indexes automatically. Each +variable uses its own available sample. The tests operate on the selected +observations without fitting a regression. + The available diagnostics are classical CD, randomized CDw, power-enhanced CDw+, and bias-corrected CD\*. CD uses pairwise-complete observations by default. CDw, CDw+, and CD\* require a balanced residual sample. Periods with no @@ -111,10 +125,24 @@ cd_test(cce, type = "all", seed = 42, ``` Use a fixed `seed` when reporting CDw or CDw+ because their Rademacher weights -are random. The tests use different corrections and should be interpreted +are random. `reps = 1` preserves the single-draw behavior. For multiple draws, +equation (33) of Juodis and Reese combines the CDw statistics as their sum +divided by `sqrt(reps)`; CDw+ adds its screening term once to that aggregate. +The authors suggest a modest number of draws, such as 30. Seeded data calls +apply the seed separately to each variable, so selection order does not affect +its result. The tests use different corrections and should be interpreted against their own assumptions; agreement among p-values is not a substitute for checking those assumptions. +CDw follows equation (30) of Juodis and Reese, using weighted covariances and +the inverse pooled variance of the demeaned input. CDw+ follows their equation +(32). CD* uses the Pesaran-Xie correction with a standardized-PCA variant: +the package standardizes units before PCA, whereas the paper uses +unstandardized observations. For fitted models, CD* removes the economic and +deterministic fitted component while retaining the common-factor component +for PCA; the other tests use full fitted residuals. Thus CD* on a fit can differ +from CD* on `residuals(fit)`. The vignette details these formulas and assumptions. + ## R model interface Fitted models support the model methods expected by downstream R tools: diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index e4c7d8a..39df9e7 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -49,13 +49,14 @@ @article{Ditzen2018 } @article{FanLiaoYao2015, - author = {Fan, Jianqing and Liao, Yuan and Yao, Jianfeng}, + author = {Fan, Jianqing and Liao, Yuan and Yao, Jiawei}, title = {Power Enhancement in High-Dimensional Cross-Section Tests}, journal = {Econometrica}, year = {2015}, volume = {83}, number = {4}, - pages = {1497--1541} + pages = {1497--1541}, + doi = {10.3982/ECTA12749} } @article{Pesaran2007, @@ -89,21 +90,33 @@ @article{Pesaran2021 } @article{JuodisReese2021, - author = {Juodis, Artūras and Reese, Søren}, - title = {The incidental parameters problem in testing for remaining cross-sectional correlation}, + author = {Juodis, Artūras and Reese, Simon}, + title = {The incidental parameters problem in testing for remaining cross-section correlation}, journal = {Journal of Business and Economic Statistics}, - year = {2021}, + year = {2022}, volume = {40}, number = {3}, - pages = {1191--1203} + pages = {1191--1203}, + doi = {10.1080/07350015.2021.1906687} } -@article{PesaranXie2021, - author = {Pesaran, M. Hashem and Xie, Yue}, - title = {A bias-corrected CD test for error cross-sectional dependence in panel models}, - journal = {Econometric Reviews}, +@techreport{PesaranXie2021, + author = {Pesaran, M. Hashem and Xie, Yimeng}, + title = {A bias-corrected CD test for error cross-sectional dependence in panel data models with latent factors}, + institution = {CESifo}, + type = {Working Paper}, + number = {9234}, year = {2021}, - volume = {41}, - number = {6}, - pages = {649--677} + note = {Revised August 2022}, + url = {https://www.ifo.de/DocDL/cesifo1_wp9234.pdf} +} + +@article{PesaranXie2026, + author = {Pesaran, M. Hashem and Xie, Yimeng}, + title = {How to detect network dependence in latent factor models? A bias-corrected CD test}, + journal = {Econometric Theory}, + year = {2026}, + pages = {1--41}, + note = {First View}, + doi = {10.1017/S0266466625100212} } diff --git a/man/cd_test.Rd b/man/cd_test.Rd index 1d2b65c..0f6204e 100644 --- a/man/cd_test.Rd +++ b/man/cd_test.Rd @@ -1,11 +1,14 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils_cd.R +% Please edit documentation in R/utils_cd.R, R/utils_cd_data.R \name{cd_test} \alias{cd_test} \alias{cd_test.default} \alias{cd_test.csdm_fit} \alias{print.cd_test} -\title{Cross-sectional dependence (CD) tests for panel residuals} +\alias{cd_test.data.frame} +\alias{cd_test.pdata.frame} +\alias{print.cd_test_list} +\title{Cross-sectional dependence (CD) tests for panel data and residuals} \usage{ cd_test(object, ...) @@ -16,6 +19,7 @@ cd_test(object, ...) seed = NULL, min_overlap = 2L, na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L, ... ) @@ -26,15 +30,47 @@ cd_test(object, ...) seed = NULL, min_overlap = 2L, na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L, ... ) \method{print}{cd_test}(x, digits = 3, ...) + +\method{cd_test}{data.frame}( + object, + ..., + id = NULL, + time = NULL, + type = c("CD", "CDw", "CDw+", "CDstar", "all"), + n_pc = 4L, + seed = NULL, + min_overlap = 2L, + na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L +) + +\method{cd_test}{pdata.frame}( + object, + ..., + id = NULL, + time = NULL, + type = c("CD", "CDw", "CDw+", "CDstar", "all"), + n_pc = 4L, + seed = NULL, + min_overlap = 2L, + na.action = c("pairwise", "drop.incomplete.times"), + reps = 1L +) + +\method{print}{cd_test_list}(x, digits = 3, ...) } \arguments{ -\item{object}{A \code{csdm_fit} model object or a numeric matrix of residuals (N x T).} +\item{object}{A \code{data.frame}, \code{pdata.frame}, \code{csdm_fit} model, +or numeric matrix with units in rows and time periods in columns (N x T).} -\item{...}{Additional arguments passed to methods.} +\item{...}{For data frames, explicitly selected numeric columns as bare names, +quoted names, or character vectors of names. Selections must be unnamed. +For matrices and fitted models, additional arguments passed to methods.} \item{type}{Which test(s) to compute: one of \code{"CD"}, \code{"CDw"}, \code{"CDw+"}, \code{"CDstar"}, or \code{"all"} (default: \code{"CD"}).} @@ -43,16 +79,27 @@ cd_test(object, ...) \item{seed}{Integer seed for weight draws. Seeded calls restore the caller's RNG state; NULL uses the current RNG stream.} -\item{min_overlap}{Minimum number of overlapping time periods required for a unit -pair to be included in CD/CDw/CDw+ (default 2).} +\item{min_overlap}{Minimum number of finite observations for retaining a unit +before time filtering (all tests), and minimum overlapping time periods +for including a pair in classical CD (default 2).} \item{na.action}{How to handle missing data: \code{"drop.incomplete.times"} removes time periods with any missing observations to create a balanced panel for CD*; \code{"pairwise"} (default) uses pairwise correlations for CD; unbalanced CDw/CDw+ requests error and CD* warns.} -\item{x}{An object of class \code{cd_test}.} +\item{reps}{Positive integer number of independent weight draws for CDw and +CDw+ (default 1). The draws are summed and divided by \code{sqrt(reps)}. +CD and CD* are unaffected. Weighted test results record the number of draws +in their \code{reps} field. For data inputs, \code{seed} is applied separately +to each variable's test, making seeded results independent of selection order.} + +\item{x}{An object of class \code{cd_test} or \code{cd_test_list}.} \item{digits}{Number of digits to print (default 3).} + +\item{id, time}{Column names (strings) for the unit and time indexes of a plain +data frame. For a \code{pdata.frame}, these are inferred from its stored +indexes and supplied values are ignored. Time labels need not be numeric.} } \value{ An object of class \code{cd_test} with fields \code{tests}, \code{type}, @@ -60,46 +107,106 @@ An object of class \code{cd_test} with fields \code{tests}, \code{type}, \code{excluded_times}, \code{kept_times}, and \code{call}. The \code{tests} list contains one or more test results, each with \code{statistic} and \code{p.value}. +Data-frame methods return a \code{cd_test_list}: a named list containing one +\code{cd_test} result per selected variable, including when only one variable +is selected. Each result also records \code{variable}, \code{units}, and +\code{excluded_unit_ids}. The list has \code{call}, \code{id}, and \code{time} +attributes and prints a combined table with each variable's sample dimensions. +For fitted models, CD* additionally records \code{input = "partial_residuals"} +and its own sample exclusions in \code{tests$CDstar}. With \code{type = "all"}, +top-level sample fields describe the full residual panel; the CD* sample +dimensions are recorded in \code{tests$CDstar$N} and \code{tests$CDstar$T}. } \description{ Computes Pesaran CD, CDw, CDw+, and CD* tests for cross-sectional dependence -in panel residuals. The implementation supports residual matrices or fitted -\code{csdm_fit} objects and provides consistent handling of unbalanced panels. +in panel variables or residuals. The implementation supports data frames, +indexed panel data frames, numeric matrices, and fitted \code{csdm_fit} objects. } \details{ +\subsection{Selecting panel variables}{ + +For a plain data frame, supply distinct unit and time column names through +\code{id} and \code{time}. For a \code{pdata.frame}, the stored indexes are +used, even when the index columns have been removed from the data. Select at +least one numeric non-index column explicitly. Each variable is tested +separately using its own available sample and the same missing-data policy. +No regression is fitted. CDw demeans the observations within each unit, as it +does for matrix inputs. Test controls and panel indexes follow \code{...} in +the data-frame methods and must be named. +} + \subsection{Notation}{ -Let \eqn{E} be the residual matrix with \eqn{N} cross-sectional units and \eqn{T} -time periods. For each unit pair \eqn{(i,j)}, let \eqn{T_{ij}} be the number of -overlapping time periods and \eqn{\rho_{ij}} the pairwise correlation. +Let \eqn{E} contain the selected observations or residuals, with \eqn{N} +cross-sectional units and \eqn{T} time periods. For each unit pair +\eqn{(i,j)}, let \eqn{T_{ij}} be the number of +overlapping time periods and \eqn{\widehat\rho_{ij}} the correlation computed +after demeaning both series over that pair's overlapping observations. } \subsection{Test statistics}{ \describe{ \item{CD (Pesaran, 2015)}{ -\deqn{CD = \sqrt{\frac{2}{N(N-1)}} \sum_{i 2\sqrt{\log(N)/T}\right\}.} -Here \eqn{\rho_{ij}} is the ordinary residual correlation, without +Equation (32) of Juodis and Reese applies the power-enhancement principle +of Fan, Liao, and Yao. Its statistic is +\deqn{CD_{W+} = CD_W + \sum_{i 2\sqrt{\log(N)/T}\right\}.} +Here \eqn{\widehat\rho_{ij}} is the ordinary residual correlation, without multiplication by \eqn{\sqrt{T}}. The nonnegative screening term is asymptotically zero under the conditions of the null hypothesis. +With \code{reps > 1}, replace \eqn{CD_W} by \eqn{\overline{CD}_W} and add +the screening term once. } -\item{CD* (Pesaran and Xie, 2021)}{ -CD is computed on residuals after removing \code{n_pc} principal components -from \eqn{E}. This provides a bias-corrected test under multifactor errors. +\item{CD* (Pesaran-Xie correction; standardized-PCA variant)}{ +Each unit is demeaned and divided by its sample standard deviation before +extracting \code{n_pc} principal components. Write \eqn{A} for this +\eqn{T\times N} standardized matrix. The estimated loading matrix +\eqn{\widehat\Gamma} is normalized so that +\eqn{\widehat\Gamma'\widehat\Gamma/N=I}, and the factors are +\eqn{\widehat F=A\widehat\Gamma/N}. For the factor-filtered residuals +\eqn{r_{it}} of \eqn{A-\widehat F\widehat\Gamma'}, define +\deqn{\widehat\sigma_i=\sqrt{T^{-1}\sum_t r_{it}^2},\quad + \widehat\varphi=N^{-1}\sum_i\widehat\gamma_i/\widehat\sigma_i,\quad + \widehat a_i=1-\widehat\sigma_i\widehat\varphi'\widehat\gamma_i.} +The implemented bias correction is +\deqn{\widehat\theta=1-N^{-1}\sum_i\widehat a_i^2,\qquad + CD^*=\frac{CD(r)+\sqrt{T/2}\,\widehat\theta}{1-\widehat\theta}.} +This uses the Pesaran-Xie correction algebra, but their PCA procedure +operates on unstandardized observations. Standardizing before PCA can +change the estimated factor space and is an implementation variant. +Setting \code{n_pc = 0} returns classical CD without factor removal. +For fitted models, the PCA input is +\eqn{\widehat v_{it}=y_{it}-x_{it}'\widehat\beta_i}, where \eqn{x_{it}} +includes all economic and deterministic design columns, including any +intercept, trend, and constructed lags. The fitted CSA contribution is +retained, following the input construction in Pesaran-Xie's regression +procedure. CD, CDw, and CDw+ use the full regression residuals instead. } } @@ -107,6 +214,16 @@ CD* requires a nondegenerate bias-correction denominator. Near-zero denominators can produce severe size distortions, including proportional loading/error-scale designs after standardization. Numerical rank checks do not establish the validity of the asymptotic approximation. + +All p-values use the two-sided standard-normal approximation +\eqn{2\Phi(-|\mathrm{statistic}|)}. The implemented tests do not estimate a +serial-correlation correction. Stationarity, time-series dependence, factor +strength, and relative panel dimensions must satisfy the relevant theory; +accepting a numeric variable does not establish these conditions. For CD, +omitting pairs or using very short overlaps can also affect calibration. +Classical CD can be biased on CCE residuals because common time parameters +have been estimated. CDw addresses this problem under Juodis-Reese's +assumptions, including \eqn{\sqrt{T}/N\to0}. } \subsection{Missing data and balance}{ @@ -138,6 +255,12 @@ cd_test(E_dep, type = "CD") # Specific test with parameters cd_test(E_indep, type = "CDstar", n_pc = 2) +# Test raw panel variables separately +panel <- expand.grid(id = 1:10, year = 1:10) +panel$x <- rnorm(nrow(panel)) +panel$y <- rnorm(nrow(panel)) +cd_test(panel, x, "y", id = "id", time = "year") + # From a fitted csdm model data(PWT_60_07, package = "csdm") df <- PWT_60_07 @@ -164,4 +287,6 @@ cd_test(fit, type = "all") \insertRef{FanLiaoYao2015}{csdm} \insertRef{PesaranXie2021}{csdm} + +\insertRef{PesaranXie2026}{csdm} } diff --git a/man/prepare_cd_input.Rd b/man/prepare_cd_input.Rd index 23475ff..733bff4 100644 --- a/man/prepare_cd_input.Rd +++ b/man/prepare_cd_input.Rd @@ -48,8 +48,9 @@ For compatibility, this function still performs: \details{ This helper is retained temporarily for compatibility. Row scaling changes the weighted covariance underlying CDw, while cross-sectional time demeaning -can mechanically induce dependence. Neither transformation is applied by -\code{cd_test()}. +can mechanically induce dependence. CDw does not standardize unit variances +or demean across units. CD* uses within-unit standardization before PCA, +as described in \code{\link[=cd_test]{cd_test()}}. \subsection{Transformation steps}{ \enumerate{ @@ -61,9 +62,9 @@ can mechanically induce dependence. Neither transformation is applied by \subsection{Why this preprocessing matters}{ -CD-type tests are sensitive to scale heterogeneity and sparse columns in -unbalanced panels. This helper creates a better-conditioned input matrix while -preserving as much usable information as possible. +These transformations change the diagnostic input. Their suitability depends +on the selected test and its assumptions; preprocessing does not establish +the validity of the test's reference distribution. } } \examples{ diff --git a/review/CD_ENHANCEMENT_VALIDATION.md b/review/CD_ENHANCEMENT_VALIDATION.md new file mode 100644 index 0000000..f5d0663 --- /dev/null +++ b/review/CD_ENHANCEMENT_VALIDATION.md @@ -0,0 +1,65 @@ +# CD enhancement validation + +Date: 2026-10-07. Branch: `feature/cd_enhancement`. +Baseline: `e562783` (`master`, csdm 2.0.0). +Environment: Windows, R 4.5.3, roxygen2 8.1.0. + +## Confirmed interface + +- Select numeric variables explicitly through `...`, using bare names, quoted + names, or character vectors. Data methods require named controls after `...`. +- Plain data frames require distinct `id` and `time` column names. Indexed panel + data frames use their stored indexes, including with `drop.index = TRUE`. +- Return one named `cd_test` result per variable in a `cd_test_list`, also for a + single selection. Print a combined table of tests and sample dimensions. +- Default `reps = 1L`. For multiple draws, CDw uses the sum of its statistics + divided by `sqrt(reps)`, following Juodis and Reese's equation (33). CDw+ adds + the screening term once to that aggregate. Seeded calls restore RNG state. +- Leave the package version at 2.0.0; describe the changes in the development + section of NEWS until the release version is chosen. + +## Checks + +- The full source testthat suite passed. Additional sample-exclusion checks + passed after the final metadata change. +- Independent cross-product references verify equation (33), the single-draw + calculation, normal p-values, and the CDw+ screening term. +- Data tests cover selection syntax, shuffled rows, structural gaps, independent + missing samples, invalid keys, labeled/date time indexes, removed pdata.frame + index columns, result printing, and RNG behavior. +- `R CMD build --no-manual .` passed, including vignette creation. +- `R CMD check --no-manual --output=.dev-check/cd-enhancement + csdm_2.0.0.tar.gz` finished with **Status: OK**, with zero errors, warnings, + and notes. All 465 installed-test assertions passed. Spelling, examples, S3 consistency, documentation, + isolated namespace loading, and vignette rebuilding all passed. + +The full package check used the normal execution context because the restricted +vignette subprocess encountered Windows Application Control DLL-loading blocks. +No Windows security settings or dependency installations were changed. + +For this host, disable renv startup with `R_PROFILE_USER=NUL` and +`R_ENVIRON_USER=NUL`, set `LANG=C` and `LC_ALL=C`, include `.dev-library` in +`R_LIBS`, and set `R_USER_CONFIG_DIR` to `.dev-check/config`. Set +`RSTUDIO_PANDOC` to the existing directory +`C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools`. +Build/check logs are retained locally under `.dev-check`. + +## Statistical smoke checks + +Run `Rscript --vanilla review/validate-cd-reps.R` from the repository root. +The script uses seed 10072026 and 500 simulations for each combination of +panel dimensions, null/alternative, and weight draw count. + +| N | T | Draws | CDw null rejection | CDw+ null rejection | +|---|---|---|---|---| +| 50 | 100 | 1 | 6.0% | 6.2% | +| 50 | 100 | 30 | 5.4% | 5.6% | +| 100 | 50 | 1 | 5.6% | 5.6% | +| 100 | 50 | 30 | 4.4% | 4.6% | + +All null designs use independent Gaussian observations with heterogeneous unit +scales. Under an added common Gaussian shock, CDw+ rejected in all 500 +simulations for each setting. Exact binomial confidence intervals and the CDw +alternative results are retained in `cd-reps-calibration.csv`. +The script passed its predeclared broad size and CDw+ power bounds. These checks +cover selected designs and do not establish calibration for every data process. diff --git a/review/CD_PAPER_AUDIT.md b/review/CD_PAPER_AUDIT.md new file mode 100644 index 0000000..79c20ac --- /dev/null +++ b/review/CD_PAPER_AUDIT.md @@ -0,0 +1,127 @@ +# CD implementation and documentation audit + +Date: 2026-10-07. Branch: `feature/cd_enhancement`. +Reviewed baseline: `aee1f29`. Package version: 2.0.0, development changes. + +## Sources and equation mapping + +| Test | Primary reference | Finding | +|---|---|---| +| Classical CD | [Pesaran's original diagnostic paper](https://docs.iza.org/dp1240.pdf), section 9, equation (31); [weak-dependence paper](https://docs.iza.org/dp6432.pdf) | The overlap-specific correlations, square-root overlap factors, and retained-unit normalization agree. | +| CDw | [Juodis-Reese, final published version](https://pure.uva.nl/ws/files/116332884/The_Incidental_Parameters_Problem_in_Testing_for_Remaining_Cross_Section_Correlation.pdf), p. 1197, equation (30) | The weighted covariance numerator and inverse pooled variance agree. | +| CDw+ | Same paper, equation (32) | The absolute-correlation screen and threshold agree; Fan-Liao-Yao provide the general enhancement principle. | +| Repeated CDw | Same paper, p. 1198, equation (33) | Sum independent draws and divide by the square root of the draw count. | +| CD* | [Pesaran-Xie working paper](https://www.ifo.de/DocDL/cesifo1_wp9234.pdf); [current author manuscript](https://arxiv.org/pdf/2109.00408), equations (12), (28)-(31), (38)-(40), Remark 6 | Correction algebra agrees. Standardization before PCA is a package variant. Fitted input needed correction. | + +The published Juodis-Reese PDF was inspected visually. Its equation (29) is +the unscaled covariance expression; equation (30) includes the inverse pooled +variance. The previous vignette already contained that factor. Its compact pair +sum and different residual symbol were equivalent, rather than a different +statistic. The revised vignette uses explicit index bounds and identifies the +equation, residual demeaning, and variance separately. Since Rademacher weights +have square one, this variance is the mean squared demeaned input. Replacing it +with pair-specific standard deviations would implement a different statistic. + +Repeated CDw+ adds one screen to the repeated CDw aggregate. This is the +explicitly agreed package extension combining equations (32) and (33). + +## CD* decisions and corrected behavior + +Two distinct choices were discussed and confirmed with the maintainer: + +1. Preserve standardization before PCA and document it as an implementation + variant. The reference PCA estimates loadings from unstandardized data; + standardization can change the factor space for heterogeneous unit variances. +2. For fitted models, start CD* with the response minus the fitted economic and + deterministic component, retaining the common-factor component. The previous + code used full CCE residuals, which is a different input from the regression + procedure described in the reference. + +The implementation reconstructs that partial residual from the stored economic +design, response, unit coefficients, and original observation map. Economic +lags, intercepts, and trends are subtracted; CSA contributions are retained. +Only estimated observations are used. CD, CDw, and CDw+ continue to test full +fitted residuals. In `type = "all"`, CD* is computed separately, with its own +sample dimensions and exclusions recorded under `tests$CDstar`. Passing the +full residual matrix explicitly to CD* retains its literal matrix-input meaning. + +The loading normalization is Gamma'Gamma/N = I. Factor-filtered residual +scales are unit-specific root mean squares. The denominator is the mean of +squared loading/scale adjustment terms, and the numerator adds the associated +square-root-time bias correction. The documentation now gives the complete +calculation. Zero requested PCs reduce to classical CD on the selected input. + +The partial-residual construction aligns with the regression input in the +reference; it does not make the retained standardized-PCA variant an exact +reproduction of the paper, nor establish the paper's assumptions for every +MG, dynamic CCE, or CS-ARDL specification. + +## Interpretation and references + +All p-values use two-sided standard-normal approximations. There is no +implemented serial-correlation adjustment. Stationarity, factor strength, +residual properties, and panel asymptotics must be assessed for the application. +Ordinary CD can be biased after fitting common time parameters. Matrix/data +acceptance is an interface feature, not evidence that a particular series meets +the asymptotic conditions. For classical CD, the minimum overlap of two is +computational; skipped pairs or very short overlaps can affect calibration. + +Corrected references: + +- Pesaran-Xie's coauthor is Yimeng Xie. The 2021 source is a working paper, + rather than the previously listed Econometric Reviews article. Retained that + working-paper reference and added the [2026 Econometric Theory publication](https://doi.org/10.1017/S0266466625100212). +- Fan-Liao-Yao's third author is Jiawei Yao, as shown in the + [published paper hosted by Fan](https://fan.princeton.edu/document/591). +- Added the verified Juodis-Reese and Fan-Liao-Yao DOIs. + +## Verification + +- Literal pair-overlap and balanced-panel references for classical CD. +- Literal triple-sum reference for CDw equation (30), including heterogeneous + scales that distinguish it from a weighted-correlation statistic. +- Existing independent cross-product references for equation (33), normal + p-values, the one-draw result, and adding the screen once. +- Independent eigenvector/loading normalization for two-factor CD*, applied to + the agreed standardized input; verified invariance to positive unit scaling. +- Fitted partial residuals checked against direct response-minus-economic-term + calculations and independent augmented regressions. Covered static CCE, + shuffled data, missing values, sample subsets, dynamic terms, CS-ARDL, trends, + MG, and numeric factor time labels in `pdata.frame` fits; other fitted test + statistics retain their full-residual references. + +These algebra checks verify computations; they are not universal size/power +guarantees for the retained PCA variant or arbitrary data processes. + +Run `Rscript --vanilla review/validate-cdstar-fit.R` for the additional fitted +null simulation. It uses seed 10072027, N=50, T=100, 500 replications, one +Gaussian common factor, heterogeneous loadings/scales, a factor-dependent +regressor, and independent Gaussian idiosyncratic errors. The correctly +specified CCE fit is followed by one-PC CD* with standardized PCA. + +| CD* input | Rejection at 5% | Exact 95% binomial interval | +|---|---|---| +| Corrected partial residuals | 5.8% | 3.92%-8.22% | +| Previous full CCE residuals | 100% | 99.26%-100% | + +The partial-residual result passed the predeclared broad 12% smoke bound. +The full-residual result is a comparison, not an acceptance case. This single +design supports the input correction but does not establish general calibration +of the standardized-PCA variant. Results are in `cdstar-fit-calibration.csv`. + +Final validation on Windows 11 with R 4.5.3: + +- Source test suite passed, followed by a successful source-package build + including the vignette. +- `R CMD check --no-manual` returned `Status: OK`: no check errors, warnings, + or notes. Installed-package tests reported 502 passes, zero failures, + zero warnings, and zero skips. +- Examples, spelling, S3 consistency, documentation, and rebuilding the + vignette all passed. The generated HTML contains the corrected equation + references, standardized-PCA qualification, partial-residual explanation, + and publication DOIs. +- The dependency check printed a temporary CRAN index download failure but + completed successfully using installed dependencies. This local check does + not verify hosted GitHub Actions startup. + +Local check logs are under `.dev-check/cd-paper-audit/` (ignored by Git). diff --git a/review/cd-reps-calibration.csv b/review/cd-reps-calibration.csv new file mode 100644 index 0000000..1aff98f --- /dev/null +++ b/review/cd-reps-calibration.csv @@ -0,0 +1,17 @@ +"N","T","alternative","reps","method","simulations","rejection","lower","upper" +50,100,FALSE,1,"CDw",500,0.06,0.0408441817627769,0.0845485537076002 +50,100,FALSE,1,"CDw_plus",500,0.062,0.0425110612963855,0.0868525974791013 +50,100,FALSE,30,"CDw",500,0.054,0.0358842597219337,0.0775974684949372 +50,100,FALSE,30,"CDw_plus",500,0.056,0.0375304444750222,0.0799213006043601 +50,100,TRUE,1,"CDw",500,0.644,0.600276004517336,0.686009857806374 +50,100,TRUE,1,"CDw_plus",500,1,0.992649389948092,1 +50,100,TRUE,30,"CDw",500,0.698,0.655670036677752,0.7379700225669 +50,100,TRUE,30,"CDw_plus",500,1,0.992649389948092,1 +100,50,FALSE,1,"CDw",500,0.056,0.0375304444750222,0.0799213006043601 +100,50,FALSE,1,"CDw_plus",500,0.056,0.0375304444750222,0.0799213006043601 +100,50,FALSE,30,"CDw",500,0.044,0.0277767672855717,0.0658610003235949 +100,50,FALSE,30,"CDw_plus",500,0.046,0.0293801753924028,0.0682254278812379 +100,50,TRUE,1,"CDw",500,0.496,0.451310742261143,0.540736812721721 +100,50,TRUE,1,"CDw_plus",500,1,0.992649389948092,1 +100,50,TRUE,30,"CDw",500,0.562,0.517253482503963,0.606009252357985 +100,50,TRUE,30,"CDw_plus",500,1,0.992649389948092,1 diff --git a/review/cdstar-fit-calibration.csv b/review/cdstar-fit-calibration.csv new file mode 100644 index 0000000..dbe8a7d --- /dev/null +++ b/review/cdstar-fit-calibration.csv @@ -0,0 +1,3 @@ +"input","seed","simulations","N","T","n_pc","rejection","lower","upper" +"partial",10072027,500,50,100,1,0.058,0.0391838687971878,0.0822382208938354 +"full",10072027,500,50,100,1,1,0.992649389948092,1 diff --git a/review/validate-cd-reps.R b/review/validate-cd-reps.R new file mode 100644 index 0000000..f3e4619 --- /dev/null +++ b/review/validate-cd-reps.R @@ -0,0 +1,35 @@ +# Reproducible size/power smoke checks for one and thirty CDw weight draws. +dir.create(".dev-check/config", recursive = TRUE, showWarnings = FALSE) +Sys.setenv(R_USER_CONFIG_DIR = normalizePath(".dev-check/config")) +.libPaths(c(normalizePath(".dev-library"), .libPaths())) +pkgload::load_all(".", quiet = TRUE) +set.seed(10072026) +simulations <- 500L +designs <- data.frame(N = c(50L, 100L), T = c(100L, 50L)) +results <- list() +for (design in seq_len(nrow(designs))) { + N <- designs$N[design] + TT <- designs$T[design] + for (alternative in c(FALSE, TRUE)) for (weight_reps in c(1L, 30L)) { + p <- replicate(simulations, { + E <- matrix(rnorm(N * TT), N) * seq(.5, 1.5, length.out = N) + if (alternative) E <- E + rep(rnorm(TT), each = N) + tests <- cd_test(E, type = "CDw+", reps = weight_reps)$tests + c(CDw = tests$CDw$p.value, CDw_plus = tests$CDw_plus$p.value) + }) + for (method in rownames(p)) { + rejects <- sum(p[method, ] < .05) + interval <- binom.test(rejects, simulations)$conf.int + results[[length(results) + 1L]] <- data.frame( + N, T = TT, alternative, reps = weight_reps, method, simulations, + rejection = rejects / simulations, lower = interval[1L], upper = interval[2L] + ) + } + } +} +out <- do.call(rbind, results) +print(out, row.names = FALSE) +write.csv(out, "review/cd-reps-calibration.csv", row.names = FALSE) +# Broad predeclared bounds detect severe failures in these selected designs. +stopifnot(all(out$rejection[!out$alternative] < .12)) +stopifnot(all(out$rejection[out$alternative & out$method == "CDw_plus"] > .8)) diff --git a/review/validate-cdstar-fit.R b/review/validate-cdstar-fit.R new file mode 100644 index 0000000..88db865 --- /dev/null +++ b/review/validate-cdstar-fit.R @@ -0,0 +1,32 @@ +# Selected null calibration for CD* on fitted partial residuals. +dir.create(".dev-check/config", recursive = TRUE, showWarnings = FALSE) +Sys.setenv(R_USER_CONFIG_DIR = normalizePath(".dev-check/config")) +.libPaths(c(normalizePath(".dev-library"), .libPaths())) +pkgload::load_all(".", quiet = TRUE) +set.seed(10072027) +simulations <- 500L +N <- 50L; TT <- 100L +d <- expand.grid(time = seq_len(TT), id = seq_len(N)) +loading <- seq(-1, 2, length.out = N) +sigma <- seq(.5, 1.5, length.out = N) +p <- replicate(simulations, { + factor <- rnorm(TT) + d$x <- .4 * factor[d$time] + rnorm(nrow(d)) + d$y <- 1 + .6 * d$x + loading[d$id] * factor[d$time] + + sigma[d$id] * rnorm(nrow(d)) + fit <- csdm(y ~ x, d, "id", "time", model = "cce") + c(partial = cd_test(fit, type = "CDstar", n_pc = 1L)$tests$CDstar$p.value, + full = cd_test(residuals(fit), type = "CDstar", n_pc = 1L)$tests$CDstar$p.value) +}) +results <- lapply(rownames(p), function(input) { + rejects <- sum(p[input, ] < .05) + interval <- binom.test(rejects, simulations)$conf.int + data.frame(input, seed = 10072027, simulations, N, T = TT, n_pc = 1L, + rejection = rejects / simulations, lower = interval[1], upper = interval[2]) +}) +out <- do.call(rbind, results) +print(out, row.names = FALSE) +write.csv(out, "review/cdstar-fit-calibration.csv", row.names = FALSE) +# Broad predeclared smoke bound for the corrected input in this selected design. +# The full-residual variant is retained as a comparison, not an acceptance case. +stopifnot(all(is.finite(p)), out$rejection[out$input == "partial"] < .12) diff --git a/tests/testthat/test-cd-data.R b/tests/testthat/test-cd-data.R new file mode 100644 index 0000000..6a259ab --- /dev/null +++ b/tests/testthat/test-cd-data.R @@ -0,0 +1,134 @@ +test_that("data-frame variables match independently shaped matrices", { + d <- panel_fixture() + E <- t(matrix(d$x, nrow = 60L)) + dimnames(E) <- list(as.character(1:8), as.character(1:60)) + expected <- cd_test(E, type = "all", n_pc = 2L, seed = 19) + shuffled <- d[sample(nrow(d)), ] + actual <- cd_test(shuffled, x, "y", id = "id", time = "time", + type = "all", n_pc = 2L, seed = 19) + expect_s3_class(actual, "cd_test_list") + expect_named(actual, c("x", "y")) + expect_s3_class(actual$x, "cd_test") + expect_equal(actual$x$tests, expected$tests) + expect_identical(actual$x$units, rownames(E)) + expect_identical(actual$x$kept_times, colnames(E)) + expect_identical(actual$x$variable, "x") + expect_identical(attr(actual, "id"), "id") + expect_identical(attr(actual, "time"), "time") + expect_true(any(grepl("shuffled", deparse(attr(actual, "call"))))) +}) + +test_that("variable selection supports bare names, strings, and caller vectors", { + d <- panel_fixture() + vars <- c("x", "y") + expected <- cd_test(d, x, y, id = "id", time = "time") + expect_equal(cd_test(d, "x", "y", id = "id", time = "time")$x$tests, + expected$x$tests) + expect_equal(cd_test(d, vars, id = "id", time = "time")$y$tests, + expected$y$tests) + expect_named(cd_test(d, c("x", "y"), id = "id", time = "time"), vars) + # Data column names take precedence over objects in the caller. + x <- "z" + expect_named(cd_test(d, x, id = "id", time = "time"), "x") + d$type <- d$x + expect_named(cd_test(d, type, id = "id", time = "time"), "type") +}) + +test_that("data inputs retain variable-specific samples and structural gaps", { + d <- panel_fixture() + d$x[d$id == 1L] <- 1 + d$x[d$time == 1L] <- NA_real_ + d <- d[!(d$id == 2L & d$time == 2L), ] + actual <- cd_test(d, x, y, id = "id", time = "time") + expect_equal(actual$x$N, 7L) + expect_equal(actual$x$T, 59L) + expect_identical(actual$x$excluded_unit_ids, "1") + expect_identical(actual$x$excluded_times, "1") + expect_equal(actual$y$N, 8L) + expect_equal(actual$y$T, 60L) + E <- matrix(NA_real_, 8L, 60L) + E[cbind(d$id, d$time)] <- d$y + expect_equal(actual$y$tests, cd_test(E)$tests) + expect_error(cd_test(d, y, id = "id", time = "time", type = "CDw"), "balanced") + expect_message(balanced <- cd_test(d, y, id = "id", time = "time", type = "CDw", + seed = 7, na.action = "drop.incomplete.times"), + "Dropped 1") + expect_equal(balanced$y$T, 59L) + expect_identical(balanced$y$excluded_times, "2") +}) + +test_that("time indexes can be labels or dates without an estimation grid", { + d <- panel_fixture(periods = 6L) + expected <- cd_test(d, x, id = "id", time = "time")$x$tests + d$time <- paste0("period-", d$time) + expect_equal(cd_test(d, x, id = "id", time = "time")$x$tests, expected) + d$time <- as.Date("2020-01-01") + rep(c(0, 2, 7, 10, 21, 40), 8L) + expect_equal(cd_test(d, x, id = "id", time = "time")$x$tests, expected) +}) + +test_that("pdata.frame methods use stored indexes including removed columns", { + skip_if_not_installed("plm") + d <- panel_fixture() + expected <- cd_test(d, x, y, id = "id", time = "time", type = "CDw", seed = 8) + for (drop_index in c(FALSE, TRUE)) { + p <- plm::pdata.frame(d, index = c("id", "time"), drop.index = drop_index) + actual <- cd_test(p, x, "y", type = "CDw", seed = 8) + expect_equal(actual$x$tests, expected$x$tests) + expect_equal(actual$y$tests, expected$y$tests) + expect_identical(actual$x$kept_times, expected$x$kept_times) + expect_identical(attr(actual, "id"), "id") + expect_equal(cd_test(p, x, id = "ignored", time = "ignored")$x$tests, + cd_test(d, x, id = "id", time = "time")$x$tests) + } +}) + +test_that("invalid selections and panel keys fail explicitly", { + d <- panel_fixture() + expect_error(cd_test(d, id = "id", time = "time"), "explicitly") + expect_error(cd_test(d, x), "distinct columns") + expect_error(cd_test(d, x, id = "id", time = "id"), "distinct columns") + expect_error(cd_test(d, "unknown", id = "id", time = "time"), "Unknown") + expect_error(cd_test(d, x, x, id = "id", time = "time"), "only once") + expect_error(cd_test(d, x, typo = "CD", id = "id", time = "time"), "unnamed") + expect_error(cd_test(d, log(x), id = "id", time = "time"), "column names") + expect_error(cd_test(d, id, id = "id", time = "time"), "indexes") + d$label <- "a" + expect_error(cd_test(d, label, id = "id", time = "time"), "numeric vectors") + expect_error(cd_test(rbind(d, d[1L, ]), x, id = "id", time = "time"), "Duplicate") + d$id[1L] <- NA + expect_error(cd_test(d, x, id = "id", time = "time"), "keys") + d <- panel_fixture() + d$time[1L] <- Inf + expect_error(cd_test(d, x, id = "id", time = "time"), "keys") + d <- panel_fixture() + d$x[1L] <- Inf + expect_error(cd_test(d, x, id = "id", time = "time"), "infinite") + expect_error(cd_test(d[FALSE, ], x, id = "id", time = "time"), "nonempty") +}) + +test_that("combined printing shows only the requested tests and each sample", { + d <- panel_fixture() + d$x[d$id == 1L] <- 0 + one <- cd_test(d, x, id = "id", time = "time") + expect_named(one, "x") + output <- capture.output(returned <- print(one)) + expect_identical(returned, one) + expect_true(any(grepl("variable", output))) + expect_true(any(grepl("7", output))) + enhanced <- cd_test(d, x, y, id = "id", time = "time", type = "CDw+", seed = 4) + output <- capture.output(print(enhanced)) + expect_equal(sum(grepl("CDw+", output, fixed = TRUE)), 2L) + expect_false(any(grepl("CDstar", output, fixed = TRUE))) + all <- cd_test(d, x, y, id = "id", time = "time", type = "all", seed = 4) + expect_equal(sum(grepl("CDstar", capture.output(print(all)), fixed = TRUE)), 2L) +}) + +test_that("seeded data-frame diagnostics preserve caller RNG and column order", { + d <- panel_fixture() + before <- .Random.seed + a <- cd_test(d, x, y, id = "id", time = "time", type = "CDw", seed = 19) + expect_identical(.Random.seed, before) + b <- cd_test(d, y, x, id = "id", time = "time", type = "CDw", seed = 19) + expect_equal(a$x$tests, b$x$tests) + expect_equal(a$y$tests, b$y$tests) +}) diff --git a/tests/testthat/test-cd-paper-reference.R b/tests/testthat/test-cd-paper-reference.R new file mode 100644 index 0000000..2d0b5ad --- /dev/null +++ b/tests/testthat/test-cd-paper-reference.R @@ -0,0 +1,44 @@ +test_that("unbalanced CD uses overlap means and Pesaran's unit normalization", { + E <- rbind(c(1, 4, 2, 7, NA, 5, 8, 3), + c(NA, 3, 6, 2, 9, 4, NA, 7), + c(2, NA, 1, 8, 3, NA, 6, 4), + c(3, 5, NA, 1, 8, 7, 2, NA)) + pairs <- combn(nrow(E), 2L) + terms <- apply(pairs, 2L, function(pair) { + observed <- is.finite(E[pair[1L], ]) & is.finite(E[pair[2L], ]) + x <- E[pair[1L], observed]; y <- E[pair[2L], observed] + x <- x - mean(x); y <- y - mean(y) + sqrt(sum(observed)) * sum(x * y) / sqrt(sum(x^2) * sum(y^2)) + }) + expected <- sqrt(2 / (nrow(E) * (nrow(E) - 1))) * sum(terms) + actual <- cd_test(E)$tests$CD + expect_equal(actual$statistic, expected) + expect_equal(actual$p.value, 2 * pnorm(abs(expected), lower.tail = FALSE)) + expect_equal(actual$pairs_used, ncol(pairs)) + balanced <- E + balanced[is.na(balanced)] <- seq_len(sum(is.na(balanced))) + expected <- sqrt(2 * ncol(balanced) / (nrow(E) * (nrow(E) - 1))) * + sum(cor(t(balanced))[upper.tri(cor(t(balanced)))]) + expect_equal(cd_test(balanced)$tests$CD$statistic, expected) +}) + +test_that("CDw uses equation 30 pooled variance rather than pair-specific correlations", { + E <- rbind(c(1, 3, 2, 5, 4, 8), c(4, 1, 8, 3, 12, 6), + c(20, 8, 4, 16, 28, 12), c(6, 12, 3, 9, 18, 15)) + N <- nrow(E); TT <- ncol(E) + U <- E - rowMeans(E) + set.seed(42) + w <- sample(c(-1, 1), N, replace = TRUE) + pooled <- sum(w^2 * U^2) / (N * TT) + numerator <- 0 + for (t in seq_len(TT)) for (i in 2:N) for (j in seq_len(i - 1L)) { + numerator <- numerator + w[i] * U[i, t] * w[j] * U[j, t] + } + expected <- pooled^(-1) * sqrt(2 / (TT * N * (N - 1))) * numerator + actual <- cd_test(E, type = "CDw", seed = 42)$tests$CDw$statistic + expect_equal(actual, expected) + correlations <- cor(t(E)) * outer(w, w) + correlation_variant <- sqrt(2 * TT / (N * (N - 1))) * + sum(correlations[upper.tri(correlations)]) + expect_gt(abs(actual - correlation_variant), .1) +}) diff --git a/tests/testthat/test-cdstar-reference.R b/tests/testthat/test-cdstar-reference.R index b162ec9..0be4cee 100644 --- a/tests/testthat/test-cdstar-reference.R +++ b/tests/testthat/test-cdstar-reference.R @@ -15,3 +15,117 @@ test_that("CDstar preserves unit-specific scales", { expect_error(cd_test(E, type = "CDstar", n_pc = 20), "n_pc") expect_error(cd_test(E, type = "CDstar", n_pc = 1.5), "integers") }) + +test_that("standardized CDstar agrees with the loading-normalized eigenvector reference", { + set.seed(485) + N <- 18L; TT <- 90L; r <- 2L + E <- matrix(rnorm(N * TT), N) * seq(.3, 2.5, length.out = N) + + outer(seq(-1, 2, length.out = N), rnorm(TT)) + + outer(rep(c(-.7, .8, 1.4), length.out = N), rnorm(TT)) + A <- scale(t(E)) + # Pesaran-Xie loading normalization, applied to the chosen standardized input. + Q <- eigen(crossprod(A), symmetric = TRUE)$vectors[, seq_len(r), drop = FALSE] + Gamma <- sqrt(N) * Q + factors <- A %*% Gamma / N + U <- A - factors %*% t(Gamma) + sigma <- sqrt(colMeans(U^2)) + phi <- colMeans(Gamma / sigma) + a <- 1 - sigma * as.numeric(Gamma %*% phi) + q <- mean(a^2) + epsilon <- sweep(U, 2L, sigma, "/") + rho <- crossprod(epsilon) / TT + cd <- sqrt(2 * TT / (N * (N - 1))) * sum(rho[upper.tri(rho)]) + expected <- (cd + sqrt(TT / 2) * (1 - q)) / q + actual <- cd_test(E, type = "CDstar", n_pc = r)$tests$CDstar + expect_equal(crossprod(Gamma) / N, diag(r), ignore_attr = TRUE) + expect_equal(actual$statistic, expected, tolerance = 1e-10) + expect_equal(actual$p.value, 2 * pnorm(abs(expected), lower.tail = FALSE)) + expect_equal(cd_test(E * seq(.5, 3, length.out = N), type = "CDstar", n_pc = r)$tests, + cd_test(E, type = "CDstar", n_pc = r)$tests, tolerance = 1e-10) + expect_equal(cd_test(E, type = "CDstar", n_pc = 0L)$tests$CDstar$statistic, + cd_test(E)$tests$CD$statistic) +}) + +test_that("fitted CDstar retains the factor component while other tests use full residuals", { + d <- panel_fixture(n = 16L, periods = 90L) + fit <- csdm(y ~ x, d, "id", "time", model = "cce") + # Construct the paper's regression input directly from the observations. + V <- residuals(fit) + V[] <- NA_real_ + for (uid in fit$meta$included_units) { + beta <- fit$coef_i[uid, ] + rows <- which(d$id == uid) + V[uid, as.character(d$time[rows])] <- d$y[rows] - beta["(Intercept)"] - beta["x"] * d$x[rows] + } + expected <- cd_test(V, type = "CDstar", n_pc = 1L) + actual <- cd_test(fit, type = "CDstar", n_pc = 1L) + expect_equal(actual$tests$CDstar$statistic, expected$tests$CDstar$statistic) + expect_equal(actual$tests$CDstar$p.value, expected$tests$CDstar$p.value) + expect_identical(actual$tests$CDstar$input, "partial_residuals") + expect_false(isTRUE(all.equal(V, residuals(fit)))) + all <- cd_test(fit, type = "all", n_pc = 1L, seed = 19, reps = 7L) + full <- cd_test(residuals(fit), type = "all", n_pc = 1L, seed = 19, reps = 7L) + expect_equal(all$tests[c("CD", "CDw", "CDw_plus")], + full$tests[c("CD", "CDw", "CDw_plus")]) + expect_equal(all$tests$CDstar, actual$tests$CDstar) + expect_false(isTRUE(all.equal(actual$tests$CDstar$statistic, full$tests$CDstar$statistic))) +}) + +test_that("partial CDstar input respects dynamic terms and original sample rows", { + set.seed(486) + d <- panel_fixture(n = 12L, periods = 90L) + d <- d[sample(nrow(d)), ] + d$y[d$id == 1L & d$time == 25] <- NA_real_ + fit <- csdm(y ~ x, d, "id", "time", model = "cs_ardl", + lr = csdm_lr(type = "ardl", ylags = 1L, xdlags = 1L), + subset = time > 2L) + V <- csdm:::.cdstar_fit_input(fit, residuals(fit)) + used <- fit$sample$row[fit$sample$used] + keys <- d[used, c("id", "time")] + # Add back the CSA fitted component as an independent reconstruction. + W <- residuals(fit) + nuisance <- setdiff(colnames(fit$augmented_matrix), colnames(fit$model_matrix)) + for (uid in fit$meta$included_units) { + rows <- which(keys$id == uid) + unit_fit <- lm.fit(fit$augmented_matrix[rows, , drop = FALSE], + model.response(model.frame(fit))[rows]) + beta <- unit_fit$coefficients[nuisance] + beta[is.na(beta)] <- 0 + W[uid, as.character(keys$time[rows])] <- + W[uid, as.character(keys$time[rows])] + + as.numeric(fit$augmented_matrix[rows, nuisance, drop = FALSE] %*% beta) + } + expect_equal(V, W, tolerance = 1e-10) + expect_identical(is.na(V), is.na(residuals(fit))) + expected <- suppressMessages(cd_test(V, type = "CDstar", n_pc = 1L, + na.action = "drop.incomplete.times")) + actual <- suppressMessages(cd_test(fit, type = "CDstar", n_pc = 1L, + na.action = "drop.incomplete.times")) + expect_equal(actual$tests$CDstar$statistic, expected$tests$CDstar$statistic) + expect_equal(actual$excluded_times, expected$excluded_times) + mg <- csdm(y ~ x, d, "id", "time", model = "mg") + expect_equal(csdm:::.cdstar_fit_input(mg, residuals(mg)), residuals(mg)) + trend <- csdm(y ~ x, d, "id", "time", model = "cce", trend = "unit") + V_trend <- csdm:::.cdstar_fit_input(trend, residuals(trend)) + for (uid in trend$meta$included_units) { + rows <- which(d$id == uid & is.finite(d$y)) + beta <- trend$coef_i[uid, ] + expected <- d$y[rows] - beta["(Intercept)"] - beta["x"] * d$x[rows] - + beta[".csdm_trend__"] * d$time[rows] + expect_equal(as.numeric(V_trend[uid, as.character(d$time[rows])]), expected) + } + old <- fit + old$model_matrix <- NULL + expect_error(cd_test(old, type = "CDstar", n_pc = 1L), "refit") +}) + +test_that("fitted CDstar maps pdata.frame factor time labels to the fitted numeric grid", { + skip_if_not_installed("plm") + d <- panel_fixture(n = 12L, periods = 90L) + numeric_fit <- csdm(y ~ x, d, "id", "time", model = "cce") + d$time <- sprintf("%03d", d$time) + panel <- plm::pdata.frame(d, index = c("id", "time")) + fit <- csdm(y ~ x, panel, "id", "time", model = "cce") + expect_equal(cd_test(fit, type = "CDstar", n_pc = 1L)$tests, + cd_test(numeric_fit, type = "CDstar", n_pc = 1L)$tests) +}) diff --git a/tests/testthat/test-diagnostic-sample.R b/tests/testthat/test-diagnostic-sample.R index 593dab3..3af0c99 100644 --- a/tests/testthat/test-diagnostic-sample.R +++ b/tests/testthat/test-diagnostic-sample.R @@ -45,3 +45,23 @@ test_that("partially observed periods retain the requested missing-data policy", expect_error(cd_test(E, type = "CDw", seed = 7), "balanced sample") expect_equal(cd_test(E, type = "CD")$T, 20L) }) + +test_that("complete-time selection records all excluded periods in input order", { + set.seed(27) + E <- matrix(rnorm(8 * 20), 8L, dimnames = list(NULL, 2001:2020)) + E[1L, 1L] <- NA_real_ + E[, 2L] <- NA_real_ + E[2L, 3L] <- NA_real_ + expect_message(a <- cd_test(E, type = "CDw", seed = 7, + na.action = "drop.incomplete.times"), "Dropped 2") + expect_identical(a$excluded_times, colnames(E)[1:3]) + expect_identical(a$kept_times, colnames(E)[-(1:3)]) + expect_equal(a$tests, cd_test(E[, -(1:3)], type = "CDw", seed = 7)$tests) + colnames(E) <- NULL + expect_message(b <- cd_test(E, na.action = "drop.incomplete.times"), "Dropped 2") + expect_identical(b$excluded_times, 1:3) + colnames(E) <- c(rep("same", 4L), as.character(5:20)) + expect_message(c <- cd_test(E, na.action = "drop.incomplete.times"), "Dropped 2") + expect_identical(c$excluded_times, rep("same", 3L)) + expect_identical(c$kept_times[1L], "same") +}) diff --git a/tests/testthat/test-weighted-cd.R b/tests/testthat/test-weighted-cd.R index d154603..72bfb3a 100644 --- a/tests/testthat/test-weighted-cd.R +++ b/tests/testthat/test-weighted-cd.R @@ -14,3 +14,79 @@ test_that("weighted CD matches the paper-defined cross-product calculation", { E[1, 1] <- NA expect_error(cd_test(E, type = "CDw+"), "balanced") }) + +test_that("repeated CDw matches equation 33 using independent cross-products", { + set.seed(731) + E <- matrix(rnorm(1200), 12L) * seq(.5, 2, length.out = 12L) + U <- E - rowMeans(E) + reps <- 30L + set.seed(84) + draws <- replicate(reps, { + w <- sample(c(-1, 1), nrow(E), replace = TRUE) + products <- crossprod(t(U * w)) + sqrt(2 / (ncol(E) * nrow(E) * (nrow(E) - 1))) * + sum(products[upper.tri(products)]) / mean(U^2) + }) + expected <- sum(draws) / sqrt(reps) + actual <- cd_test(E, type = "CDw", seed = 84, reps = reps)$tests$CDw + expect_equal(actual$statistic, expected, tolerance = 1e-12) + expect_equal(actual$p.value, 2 * pnorm(abs(expected), lower.tail = FALSE)) + expect_identical(actual$reps, reps) + expect_equal(cd_test(E, type = "CDw", seed = 84, reps = 1L)$tests$CDw$statistic, + draws[1L], tolerance = 1e-12) + expect_identical(cd_test(E, type = "CDw", seed = 84)$tests, + cd_test(E, type = "CDw", seed = 84, reps = 1L)$tests) +}) + +test_that("repeated CDw-plus adds its screening statistic exactly once", { + set.seed(732) + E <- matrix(rnorm(12 * 100), 12L) + rep(rnorm(100), each = 12L) + rho <- abs(cor(t(E))[upper.tri(cor(t(E)))]) + enhancement <- sum(rho[rho > 2 * sqrt(log(nrow(E)) / ncol(E))]) + expect_gt(enhancement, 0) + a <- cd_test(E, type = "CDw+", seed = 18, reps = 30L) + expected <- cd_test(E, type = "CDw", seed = 18, reps = 30L)$tests$CDw$statistic + + enhancement + expect_equal(a$tests$CDw_plus$statistic, expected) + expect_equal(a$tests$CDw_plus$p.value, 2 * pnorm(abs(expected), lower.tail = FALSE)) + expect_identical(a$tests$CDw_plus$reps, 30L) + expect_equal(a$tests$CDw_plus$enhancement, enhancement) + all <- cd_test(E, type = "all", seed = 18, reps = 30L, n_pc = 1L) + expect_equal(all$tests$CDw, a$tests$CDw) + expect_equal(all$tests$CDw_plus, a$tests$CDw_plus) + expect_equal(all$tests$CD, cd_test(E)$tests$CD) + expect_equal(all$tests$CDstar, + cd_test(E, type = "CDstar", n_pc = 1L)$tests$CDstar) +}) + +test_that("repetitions are validated and forwarded through all methods", { + d <- panel_fixture() + fit <- csdm(y ~ x, d, "id", "time") + E <- residuals(fit) + expected <- cd_test(E, type = "CDw+", reps = 7L, seed = 18)$tests + expect_equal(cd_test(fit, type = "CDw+", reps = 7L, seed = 18)$tests, expected) + d$e <- as.numeric(t(E)) + expect_equal(cd_test(d, e, id = "id", time = "time", type = "CDw+", + reps = 7L, seed = 18)$e$tests, expected) + for (invalid in list(0, -1, 1.5, NA_real_, Inf, c(1, 2), "30", TRUE, NULL)) { + expect_error(cd_test(E, type = "CDw", reps = invalid), "reps") + } + skip_if_not_installed("plm") + p <- plm::pdata.frame(d, index = c("id", "time")) + expect_equal(cd_test(p, e, type = "CDw+", reps = 7L, seed = 18)$e$tests, expected) +}) + +test_that("repeated weighted tests preserve seeded RNG and consume unseeded draws", { + set.seed(733) + E <- matrix(rnorm(600), 10L) + before <- .Random.seed + cd_test(E, type = "CDw+", reps = 7L, seed = 9) + expect_identical(.Random.seed, before) + expect_error(cd_test(E, type = "all", reps = 7L, seed = 9, n_pc = 100L), "n_pc") + expect_identical(.Random.seed, before) + cd_test(E, type = "CDw", reps = 7L) + after <- .Random.seed + .Random.seed <<- before + replicate(7L, sample(c(-1, 1), nrow(E), replace = TRUE)) + expect_identical(.Random.seed, after) +}) diff --git a/vignettes/introduction_to_csdm.Rmd b/vignettes/introduction_to_csdm.Rmd index bc907c0..c461fcc 100644 --- a/vignettes/introduction_to_csdm.Rmd +++ b/vignettes/introduction_to_csdm.Rmd @@ -273,7 +273,7 @@ modelsummary::modelsummary( ) ``` -## Residual cross-sectional dependence +## Cross-sectional dependence in residuals and variables Let $e_{it}$ denote the fitted residual and $T_{ij}$ the number of overlapping finite observations for units $i$ and $j$. The classical statistic implemented @@ -290,20 +290,42 @@ approximation. Depending on the theoretical formulation, that null is stated as cross-sectional independence or sufficiently weak dependence [@Pesaran2015; @Pesaran2021]. -For a balanced residual matrix, CDw draws one independent Rademacher weight -$w_i\in\{-1,1\}$ per unit and computes [@JuodisReese2021] +For a balanced panel, CDw draws one independent Rademacher weight +$w_i\in\{-1,1\}$ per unit, independently of the data. Write +$u_{it}=e_{it}-\bar e_i$ for the input demeaned within unit. Equation (30) of +the published Juodis-Reese paper [@JuodisReese2021, p. 1197] is $$ -CD_W= -\left(\frac{1}{NT}\sum_{i,t}w_i^2\widetilde e_{it}^2\right)^{-1} +\widehat s_W^2=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}w_i^2u_{it}^2, +\qquad +CD_W=(\widehat s_W^2)^{-1} \sqrt{\frac{2}{TN(N-1)}} -\sum_t\sum_{i