Nothing
R/testcorr.R
In testcorr: Testing Zero Correlation
Defines functions
iid.test
print.ac.test
ac.test
Documented in
ac.test
iid.test
########################################################################################################################
### UNIVARIATE STATISTICS
########################################################################################################################
#' Testing zero autocorrelation
#'
#' The function ac.test computes the test statistics for examining the null hypothesis
#' of zero autocorrelation for univariate time series given in Dalla, Giraitis and Phillips (2022).
#'
#' @usage
#' ac.test(x, max.lag, m0 = 1, alpha = 0.05, lambda = 2.576,
#' plot = TRUE, var.name = NULL, scale.font = 1)
#' @param x A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param max.lag Maximum lag at which to calculate the test statistics.
#' @param m0 Minimum lag at which to calculate the cumulative test statistics. Default is 1.
#' @param alpha Significance level for hypothesis testing used in the plots. Default is 0.05.
#' @param lambda Threshold in \eqn{\widetilde{Q}}{Q-tilde} test statistics. Default is 2.576.
#' @param plot Logical. If TRUE, 1) the sample autocorrelations with their confidence bands are plotted and 2) the cumulative test statistics with their critical values are plotted. Default is TRUE. Can be a logical vector for each of the plots 1)-2).
#' @param var.name NULL or a character string specifying the variable name. If NULL and x has name, the name of x is used. If NULL and x has no name, the string "x" is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The standard \eqn{t} and robust \eqn{\widetilde{t}}{t-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_k=0}{H[0]:\rho[k]=0} at lags \eqn{k=1,...,max.lag},
#' and the standard \eqn{LB} and robust \eqn{\widetilde{Q}}{Q-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_{m_0}=...=\rho_m=0}{H[0]:\rho[m0]=...=\rho[m]=0} at lags \eqn{m=m_0,...,max.lag}{m=m0,...,max.lag},
#' where \eqn{\rho_k}{\rho[k]} denotes the autocorrelation of \eqn{x_t}{x[t]} at lag \eqn{k}.
#' @return An object of class "ac.test", which is a list with the following components:
#' \item{lag}{The lags of the sample autocorrelations.}
#' \item{ac}{The sample autocorrelations.}
#' \item{scb}{The lower and upper limit of the confidence bands based on the standard test statistics.}
#' \item{rcb}{The lower and upper limit of the confidence bands based on the robust test statistics.}
#' \item{t}{The \eqn{t}{t} test statistics.}
#' \item{pvt}{The p-values for the \eqn{t}{t} test statistics.}
#' \item{ttilde}{The \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{pvttilde}{The p-values for the \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{lagc}{The lags of the cumulative test statistics.}
#' \item{lb}{The \eqn{LB} test statistics.}
#' \item{pvlb}{The p-values for the \eqn{LB} test statistics.}
#' \item{qtilde}{The \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{pvqtilde}{The p-values for the \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{alpha}{Significance level for hypothesis testing used in the plots.}
#' \item{varname}{The variable name used in the plots/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}. \cr
#' Giraitis, L., Li, Y. and Phillips, P. C. B. (2024). "Robust Inference on Correlation under General Heterogeneity". Journal of Econometrics, 244(1), 105691, \doi{doi:10.1016/j.jeconom.2024.105691}.
#' @examples
#' x <- rnorm(100)
#' ac.test(x, max.lag = 10)
#' @importFrom stats acf pchisq pnorm qchisq qnorm is.ts
#' @import ggplot2
#' @importFrom scales breaks_pretty
#' @importFrom knitr kable
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
ac.test <- function(x, max.lag, m0 = 1, alpha = 0.05, lambda = 2.576, plot = TRUE, var.name = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (!is.numeric(max.lag)) stop('argument "max.lag" must be numeric')
if (!is.numeric(m0)) stop('argument "m0" must be numeric')
if (!is.numeric(alpha)) stop('argument "alpha" must be numeric')
if (!is.numeric(lambda)) stop('argument "lambda" must be numeric')
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.name) & NROW(var.name) == 1 & !is.character(var.name)) stop('argument "var.name" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) != 1) stop('argument "x" must be univariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (max.lag != ceiling(max.lag) | max.lag < m0 | max.lag > (NROW(x) - 1)) stop('argument "max.lag" must be an integer value greater or equal to "m0" and less than the sample size')
if (m0 != ceiling(m0) | m0 < 1 | m0 > max.lag) stop('argument "m0" must be an integer value greater or equal to one and less or equal than "maxlag"')
if (alpha < 0 | alpha > 1) stop('argument "alpha" must be between 0 and 1')
if (lambda < 0) stop('argument "lambda" must be non-negative')
if (is.null(var.name) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop('argument "x" must have one name')
}
if (!is.null(var.name)) {
if (NROW(var.name) != NCOL(x)) stop('argument "var.name" must contain one name')
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
if (is.null(var.name)) {
if (!is.null(colnames(x))) {
my.name <- colnames(x)
}
if (is.null(colnames(x))) {
my.name <- "x"
}
}
if (!is.null(var.name)) {
my.name <- var.name
}
x <- as.numeric(x)
n <- NROW(x)
x.tilde <- x - mean(x)
seq1maxlag <- seq(1, max.lag)
seqm0maxlag <- seq(m0, max.lag)
matrix.e <- x.tilde[2:n] * x.tilde[1:(n - 1)]
if (max.lag >= 2) {
for (kk in 2:max.lag) {
matrix.e <- cbind(matrix.e, c(matrix(0, (kk - 1), 1), x.tilde[(kk + 1):n] * x.tilde[1:(n - kk)]))
}
}
matrix.esq <- matrix.e ^ 2
if (max.lag == 1) {
t.tilde <- sum(matrix.e) / sqrt(sum(matrix.esq))
}
if (max.lag >= 2) {
t.tilde <- colSums(matrix.e) / sqrt(colSums(matrix.esq))
}
big.rstar <- diag(max.lag) / 2
if (max.lag >= 2) {
for (jj in 1:(max.lag - 1)) {
for (kk in (jj + 1):max.lag) {
max.jk <- max(jj, kk)
num.r <- t(matrix.e[, jj]) %*% matrix.e[, kk]
rjkstar.ind <- num.r / sqrt(t(matrix.esq[max.jk:(n - 1), jj]) %*% matrix.esq[max.jk:(n - 1), kk])
if (abs(rjkstar.ind) > lambda) {
big.rstar[jj, kk] <- num.r / sqrt(sum(matrix.esq[max.jk:(n - 1), jj]) * sum(matrix.esq[max.jk:(n - 1), kk]))
}
}
}
}
big.rstar <- big.rstar + t(big.rstar)
q.tilde <- matrix(NA, nrow = max.lag - m0 + 1, ncol = 1)
for (mm in m0:max.lag) {
q.tilde[mm - m0 + 1] <- t(t.tilde[m0:mm]) %*% solve(big.rstar[m0:mm, m0:mm]) %*% t.tilde[m0:mm]
}
ac <- acf(x, max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
t <- sqrt(n) * ac
sq.ac <- ac ^ 2
seq.overnmink <- (n - seq1maxlag) ^ (-1)
lb <- n * (n + 2) * cumsum(sq.ac[m0:max.lag] * seq.overnmink[m0:max.lag])
z.cv <- qnorm(1 - alpha / 2)
s.cb <- z.cv / sqrt(n) * rep(1, max.lag)
r.cb <- z.cv * ac / t.tilde
lu.s.cb <- cbind(-s.cb, s.cb)
lu.r.cb <- cbind(-r.cb, r.cb)
rownames(lu.r.cb) <- rownames(lu.s.cb)
colnames(lu.r.cb) <- colnames(lu.s.cb)
pv.t <- 2 * (1 - pnorm(abs(t)))
pv.ttilde <- 2 * (1 - pnorm(abs(t.tilde)))
pv.lb <- 1 - pchisq(lb, seqm0maxlag - m0 + 1)
pv.qtilde <- 1 - pchisq(q.tilde, seqm0maxlag - m0 + 1)
if (NROW(plot) == 1) {
plotA <- plot
plotQ <- plot
}
if (NROW(plot) == 2) {
plotA <- plot[1]
plotQ <- plot[2]
}
if (plotA == TRUE & max.lag >= 2) {
myplotcorr(max.lag, seq1maxlag, ac, s.cb, r.cb, alpha, n, 1, my.name, scale.font)
}
if (plotQ == TRUE & max.lag >= m0 + 1) {
chisq.cv <- qchisq(alpha, seqm0maxlag - m0 + 1, lower.tail = FALSE)
Sys.sleep(2)
myplotstat(max.lag, seqm0maxlag, lb, q.tilde, alpha, chisq.cv, n, 1, my.name, " LB", expression(~tilde(Q)), 1, scale.font)
}
colnames(seq1maxlag) <- NULL
rownames(seq1maxlag) <- NULL
colnames(ac) <- NULL
rownames(ac) <- NULL
colnames(lu.s.cb) <- NULL
rownames(lu.s.cb) <- NULL
colnames(lu.r.cb) <- NULL
rownames(lu.r.cb) <- NULL
colnames(t) <- NULL
rownames(t) <- NULL
colnames(pv.t) <- NULL
rownames(pv.t) <- NULL
colnames(t.tilde) <- NULL
rownames(t.tilde) <- NULL
colnames(pv.ttilde) <- NULL
rownames(pv.ttilde) <- NULL
colnames(seqm0maxlag) <- NULL
rownames(seqm0maxlag) <- NULL
colnames(lb) <- NULL
rownames(lb) <- NULL
colnames(pv.lb) <- NULL
rownames(pv.lb) <- NULL
colnames(q.tilde) <- NULL
rownames(q.tilde) <- NULL
colnames(pv.qtilde) <- NULL
rownames(pv.qtilde) <- NULL
ac.test.out <- structure(list(lag = as.numeric(seq1maxlag), ac = as.numeric(ac), scb = as.matrix(lu.s.cb), rcb = as.matrix(lu.r.cb), t = as.numeric(t), pvt = as.numeric(pv.t), ttilde = as.numeric(t.tilde), pvttilde = as.numeric(pv.ttilde), lagc = as.numeric(seqm0maxlag), lb = as.numeric(lb), pvlb = as.numeric(pv.lb), qtilde = as.numeric(q.tilde), pvqtilde = as.numeric(pv.qtilde), alpha = alpha, varname = my.name), class = "ac.test")
return(ac.test.out)
}
#' @export
print.ac.test <- function(x, ...) {
options(max.print = 1000000)
results.t <- cbind(x$t, x$pvt, x$ttilde, x$pvttilde)
results.t <- format(round(results.t, 3), nsmall = 3)
m0 <- min(x$lagc)
if (m0 == 1) {
results.q <- cbind(x$lb, x$pvlb, x$qtilde, x$pvqtilde)
}
if (m0 >= 2) {
results.q <- rbind(matrix(NA, nrow = m0 - 1, ncol = 4), cbind(x$lb, x$pvlb, x$qtilde, x$pvqtilde))
}
results.q <- format(round(results.q, 3), nsmall = 3)
results.ac <- x$ac
results.ac <- format(round(results.ac, 3), nsmall = 3)
results.cb.lu <- cbind(x$scb, x$rcb)
results.cb.lu <- format(round(results.cb.lu, 3), nsmall = 3)
results.scb <- paste("(", results.cb.lu[ , 1], ",", results.cb.lu[ , 2], ")", sep = "")
results.rcb <- paste("(", results.cb.lu[ , 3], ",", results.cb.lu[ , 4], ")", sep = "")
results.tq <- cbind(x$lag, results.ac, results.scb, results.rcb, x$lag, results.t, x$lag, results.q)
results.tq <- data.frame(results.tq, row.names = NULL)
names(results.tq) <- c("Lag", "AC", paste("Stand. CB(", 100 * (1 - x$alpha), "%)", sep = ""), paste("Robust CB(", 100 * (1 - x$alpha), "%)", sep = ""), "Lag", "t", "p-value", "t-tilde", "p-value", "Lag", "LB", "p-value", "Q-tilde", "p-value")
msg <- paste("Tests for zero autocorrelation of ", x$varname, sep = "")
table.tq <- kable(results.tq, align = "r", format = "pipe")
cat("\n")
cat(msg)
print(table.tq)
cat("\n")
}
########################################################################################################################
### UNIVARIATE STATISTICS: iid
########################################################################################################################
#' Testing iid property
#'
#' The function iid.test computes the test statistics for examining the null hypothesis
#' of i.i.d. property for univariate series given in Dalla, Giraitis and Phillips (2022).
#'
#' @usage
#' iid.test(x, max.lag, m0 = 1, alpha = 0.05,
#' plot = TRUE, var.name = NULL, scale.font = 1)
#' @param x A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param max.lag Maximum lag at which to calculate the test statistics.
#' @param m0 Minimum lag at which to calculate the cumulative test statistics. Default is 1.
#' @param alpha Significance level for hypothesis testing used in the plots. Default is 0.05.
#' @param plot Logical. If TRUE, 1) the test statistics (J) and their critical values are plotted and 2) the cumulative test statistics (C) with their critical values are plotted. Default is TRUE. Can be a logical vector for each of the plots 1)-2).
#' @param var.name NULL or a character string specifying the variable name. If NULL and x has name, the name of x is used. If NULL and x has no name, the string "x" is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The \eqn{J_{x,|x|}}{J[x,|x|]} and \eqn{J_{x,x^2}}{J[x,x^2]} statistics are for testing the null hypothesis of i.i.d. at lag \eqn{k}, \eqn{k=1,...,max.lag},
#' and the \eqn{C_{x,|x|}}{C[x,|x|]} and \eqn{C_{x,x^2}}{C[x,x^2]} statistics are for testing the null hypothesis of i.i.d. at lags \eqn{m_0,...,m}, \eqn{m=m_0,...,max.lag}.
#' @return An object of class "iid.test", which is a list with the following components:
#' \item{lag}{The lags of the test statistics.}
#' \item{jab}{The \eqn{J_{x,|x|}}{J[x,|x|]} test statistics.}
#' \item{pvjab}{The p-values for the \eqn{J_{x,|x|}}{J[x,|x|]} test statistics.}
#' \item{jsq}{The \eqn{J_{x,x^2}}{J[x,x^2]} test statistics.}
#' \item{pvjsq}{The p-values for the \eqn{J_{x,x^2}}{J[x,x^2]} test statistics.}
#' \item{lagc}{The lags of the cumulative test statistics.}
#' \item{cab}{The \eqn{C_{x,|x|}}{C[x,|x|]} test statistics.}
#' \item{pvcab}{The p-values for the \eqn{C_{x,|x|}}{C[x,|x|]} test statistics.}
#' \item{csq}{The \eqn{C_{x,x^2}}{C[x,x^2]} test statistics.}
#' \item{pvcsq}{The p-values for the \eqn{C_{x,x^2}}{C[x,x^2]} test statistics.}
#' \item{alpha}{Significance level for hypothesis testing used in the plots.}
#' \item{varname}{The variable name used in the plots/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}.
#' @examples
#' x <- rnorm(100)
#' iid.test(x, max.lag = 10)
#' @importFrom stats acf pchisq pnorm qchisq qnorm is.ts
#' @import ggplot2
#' @importFrom scales breaks_pretty
#' @importFrom knitr kable
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
iid.test <- function(x, max.lag, m0 = 1, alpha = 0.05, plot = TRUE, var.name = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (!is.numeric(max.lag)) stop('argument "max.lag" must be numeric')
if (!is.numeric(m0)) stop('argument "m0" must be numeric')
if (!is.numeric(alpha)) stop('argument "alpha" must be numeric')
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.name) & NROW(var.name) == 1 & !is.character(var.name)) stop('argument "var.name" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) != 1) stop('argument "x" must be univariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (max.lag != ceiling(max.lag) | max.lag < m0 | max.lag > (NROW(x) - 1)) stop('argument "max.lag" must be an integer value greater or equal to "m0" and less than the sample size')
if (m0 != ceiling(m0) | m0 < 1 | m0 > max.lag) stop('argument "m0" must be an integer value greater or equal to one and less or equal than "maxlag"')
if (alpha < 0 | alpha > 1) stop('argument "alpha" must be between 0 and 1')
if (is.null(var.name) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop('argument "x" must have one name')
}
if (!is.null(var.name)) {
if (NROW(var.name) != NCOL(x)) stop('argument "var.name" must contain one name')
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
if (is.null(var.name)) {
if (!is.null(colnames(x))) {
my.name <- colnames(x)
}
if (is.null(colnames(x))) {
my.name <- "x"
}
}
if (!is.null(var.name)) {
my.name <- var.name
}
x <- as.numeric(x)
n <- NROW(x)
x.tilde <- x - mean(x)
seq1maxlag <- seq(1, max.lag)
seqm0maxlag <- seq(m0, max.lag)
ac <- acf(x, max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
ac.abs <- acf(abs(x.tilde), max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
ac.sq <- acf(x.tilde ^ 2, max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
j.abs <- (n ^ 2) / (n - seq1maxlag) * (ac ^ 2 + ac.abs ^ 2)
j.sq <- (n ^ 2) / (n - seq1maxlag) * (ac ^ 2 + ac.sq ^ 2)
c.abs <- cumsum(j.abs[m0:max.lag])
c.sq <- cumsum(j.sq[m0:max.lag])
pv.jabs <- 1 - pchisq(j.abs, 2)
pv.jsq <- 1 - pchisq(j.sq, 2)
pv.cabs <- 1 - pchisq(c.abs, 2 * (seqm0maxlag - m0 + 1))
pv.csq <- 1 - pchisq(c.sq, 2 * (seqm0maxlag - m0 + 1))
if (NROW(plot) == 1) {
plotJ <- plot
plotC <- plot
}
if (NROW(plot) == 2) {
plotJ <- plot[1]
plotC <- plot[2]
}
if (plotJ == TRUE & max.lag >= 2) {
chisq2.cv <- rep(qchisq(alpha, 2, lower.tail = FALSE), max.lag)
fvar.name <- substr(my.name, 1, 1)
sub.abs <- bquote(.(fvar.name) * ",|" * .(fvar.name) * "|")
sub.sq <- bquote(.(fvar.name) * "," * .(fvar.name)^2)
name.j.abs <- bquote(~J[.(sub.abs)])
name.j.sq <- bquote(~J[.(sub.sq)])
myplotstat(max.lag, seq1maxlag, j.abs, j.sq, alpha, chisq2.cv, n, 1, my.name, name.j.abs, name.j.sq, 2, scale.font)
}
if (plotC == TRUE & max.lag >= m0 + 1) {
chisq.cv <- qchisq(alpha, 2 * (seqm0maxlag - m0 + 1), lower.tail = FALSE)
fvar.name <- substr(my.name, 1, 1)
sub.abs <- bquote(.(fvar.name) * ",|" * .(fvar.name) * "|")
sub.sq <- bquote(.(fvar.name) * "," * .(fvar.name)^2)
name.c.abs <- bquote(~C[.(sub.abs)])
name.c.sq <- bquote(~C[.(sub.sq)])
Sys.sleep(2)
myplotstat(max.lag, seqm0maxlag, c.abs, c.sq, alpha, chisq.cv, n, 1, my.name, name.c.abs, name.c.sq, 3, scale.font)
}
colnames(seq1maxlag) <- NULL
rownames(seq1maxlag) <- NULL
colnames(j.abs) <- NULL
rownames(j.abs) <- NULL
colnames(pv.jabs) <- NULL
rownames(pv.jabs) <- NULL
colnames(j.sq) <- NULL
rownames(j.sq) <- NULL
colnames(pv.jsq) <- NULL
rownames(pv.jsq) <- NULL
colnames(seqm0maxlag) <- NULL
rownames(seqm0maxlag) <- NULL
colnames(c.abs) <- NULL
rownames(c.abs) <- NULL
colnames(pv.cabs) <- NULL
rownames(pv.cabs) <- NULL
colnames(c.sq) <- NULL
rownames(c.sq) <- NULL
colnames(pv.csq) <- NULL
rownames(pv.csq) <- NULL
iid.test.out <- structure(list(lag = as.numeric(seq1maxlag), jabs = as.numeric(j.abs), pvjabs = as.numeric(pv.jabs), jsq = as.numeric(j.sq), pvjsq = as.numeric(pv.jsq), lagc = as.numeric(seqm0maxlag), cabs = as.numeric(c.abs), pvcabs = as.numeric(pv.cabs), csq = as.numeric(c.sq), pvcsq = as.numeric(pv.csq), alpha = alpha, varname = my.name), class = "iid.test")
return(iid.test.out)
}
#' @export
print.iid.test <- function(x, ...) {
options(max.print = 1000000)
results.j <- cbind(x$jabs, x$pvjabs, x$jsq, x$pvjsq)
results.j <- format(round(results.j, 3), nsmall = 3)
m0 <- min(x$lagc)
if (m0 == 1) {
results.c <- cbind(x$cabs, x$pvcabs, x$csq, x$pvcsq)
}
if (m0 >= 2) {
results.c <- rbind(matrix(NA, nrow = m0 - 1, ncol = 4), cbind(x$cabs, x$pvcabs, x$csq, x$pvcsq))
}
results.c <- format(round(results.c, 3), nsmall = 3)
results.jc <- cbind(x$lag, results.j, x$lag, results.c)
results.jc <- data.frame(results.jc, row.names = NULL)
fvar.name <- substr(x$varname, 1, 1)
name.j.abs <- paste("J[", fvar.name, ",", "\u01C0", fvar.name, "\u01C0]", sep = "")
name.j.sq <- paste("J[", fvar.name, ",", fvar.name, "\u00B2]", sep = "")
name.c.abs <- paste("C[", fvar.name, ",", "\u01C0", fvar.name, "\u01C0]", sep = "")
name.c.sq <- paste("C[", fvar.name, ",", fvar.name, "\u00B2]", sep = "")
names(results.jc) <- c("Lag", name.j.abs, "p-value", name.j.sq, "p-value", "Lag", name.c.abs, "p-value", name.c.sq, "p-value")
msg <- paste("Tests for i.i.d. property of ", x$varname, sep = "")
table.jc <- kable(results.jc, align = "r", format = "pipe")
cat("\n")
cat(msg)
print(table.jc)
cat("\n")
}
########################################################################################################################
### BIVARIATE STATISTICS
########################################################################################################################
#' Testing zero cross-correlation
#'
#' The function cc.test computes the test statistics for examining the null hypothesis
#' of zero cross-correlation for bivariate time series given in Dalla, Giraitis and Phillips (2022).
#'
#' @usage
#' cc.test(x, y, max.lag, m0 = 0, alpha = 0.05, lambda = 2.576,
#' plot = TRUE, var.names = NULL, scale.font = 1)
#' @param x A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param y A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param max.lag Maximum lag at which to calculate the test statistics.
#' @param m0 Minimum lag at which to calculate the cumulative test statistics. Default is 0.
#' @param alpha Significance level for hypothesis testing used in the plots. Default is 0.05.
#' @param lambda Threshold in \eqn{\widetilde{Q}}{Q-tilde} test statistics. Default is 2.576.
#' @param plot Logical. If TRUE, 1) the sample cross-correlations with their confidence bands are plotted and 2) the cumulative test statistics with their critical values are plotted. Default is TRUE. Can be a logical vector for each of the plots 1)-2).
#' @param var.names NULL or a character string specifying the variable names. If NULL and x,y have names, the names of x,y are used. If NULL and x,y have no names, the string c("x","y") is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The standard \eqn{t} and robust \eqn{\widetilde{t}}{t-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_k=0}{H[0]:\rho[k]=0} at lags \eqn{k=-max.lag,...,-1,0,1,max.lag},
#' and the standard \eqn{HB} and robust \eqn{\widetilde{Q}}{Q-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_{m_0}=...=\rho_m=0}{H[0]:\rho[m0]=...=\rho[m]=0} at lags \eqn{m=-max.lag,...,-1,0,1,max.lag},
#' where \eqn{\rho_k}{\rho[k]} denotes the cross-correlation of \eqn{x_t}{x[t]} and \eqn{y_{t-k}}{y[t-k]} at lag \eqn{k}.
#' @return An object of class "cc.test", which is a list with the following components:
#' \item{lag}{The lags of the sample cross-correlations.}
#' \item{cc}{The sample cross-correlations.}
#' \item{scb}{The lower and upper limit of the confidence bands based on the standard test statistics.}
#' \item{rcb}{The lower and upper limit of the confidence bands based on the robust test statistics.}
#' \item{t}{The \eqn{t}{t} test statistics.}
#' \item{pvt}{The p-values for the \eqn{t}{t} test statistics.}
#' \item{ttilde}{The \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{pvtttilde}{The p-values for the \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{lagc}{The lags of the cumulative test statistics.}
#' \item{hb}{The \eqn{HB} test statistics.}
#' \item{pvhb}{The p-values for the \eqn{HB} test statistics.}
#' \item{qtilde}{The \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{pvqtilde}{The p-values for the \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{alpha}{Significance level for hypothesis testing used in the plots.}
#' \item{varnames}{The variable names used in the plots/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}. \cr
#' Giraitis, L., Li, Y. and Phillips, P. C. B. (2024). "Robust Inference on Correlation under General Heterogeneity". Journal of Econometrics, 244(1), 105691, \doi{doi:10.1016/j.jeconom.2024.105691}.
#' @examples
#' x <- rnorm(100)
#' y <- rnorm(100)
#' cc.test(x, y, max.lag = 10)
#' @importFrom stats acf ccf pchisq pnorm qchisq qnorm is.ts
#' @import ggplot2
#' @importFrom scales breaks_pretty
#' @importFrom knitr kable
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
cc.test <- function(x, y, max.lag, m0 = 0, alpha = 0.05, lambda = 2.576, plot = TRUE, var.names = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.null(y)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (is.ts(y) | is.xts(y) | is.zoo(y) | is.numeric(y) | is.data.frame(y)) {}
else{
stop('argument "y" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (!is.numeric(max.lag)) stop('argument "max.lag" must be numeric')
if (!is.numeric(m0)) stop('argument "m0" must be numeric')
if (!is.numeric(alpha)) stop('argument "alpha" must be numeric')
if (!is.numeric(lambda)) stop('argument "lambda" must be numeric')
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.names) & NROW(var.names) == 2 & (!is.character(var.names[1]) | !is.character(var.names[2]))) stop('argument "var.names" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) != 1) stop('argument "x" must be univariate')
if (NCOL(y) != 1) stop('argument "y" must be univariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (any(is.na(y))) stop('argument "y" must not contain missing values')
if (NCOL(x) != NCOL(y)) stop('arguments "x" and "y" must have the same length')
if (max.lag != ceiling(max.lag) | max.lag < m0 | max.lag > (NROW(x) - 1)) stop('argument "max.lag" must be an integer value greater or equal to "m0" and less than the sample size')
if (m0 != ceiling(m0) | m0 < 0 | m0 > max.lag) stop('argument "m0" must be an integer value greater or equal to zero and less or equal than "maxlag"')
if (alpha < 0 | alpha > 1) stop('argument "alpha" must be between 0 and 1')
if (lambda < 0) stop('argument "lambda" must be non-negative')
if (is.null(var.names) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop('argument "x" must have one name')
}
if (is.null(var.names) & !is.null(colnames(y))) {
if (NROW(colnames(y)) != NCOL(x)) stop('argument "y" must have one name')
}
if (!is.null(var.names)) {
if (NROW(var.names) != 2) stop('argument "var.names" must contain two names')
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
if (is.null(var.names)) {
if (!is.null(colnames(x))) {
my.names <- colnames(x)
}
if (is.null(colnames(x))) {
my.names <- c("x","y")
}
}
if (!is.null(var.names)) {
my.names <- var.names
}
x <- as.numeric(x)
y <- as.numeric(y)
n <- NROW(x)
results.tq.xyk <- cc.test.t.tmink(x, y, max.lag, m0, alpha, lambda)
results.tq.yxk <- cc.test.t.tmink(y, x, max.lag, m0, alpha, lambda)
results.tq <- rbind(results.tq.yxk[(max.lag + 1):2, ], results.tq.xyk)
seqmaxlagmaxlag <- seq(-max.lag, max.lag)
if (m0 == 0) {
seqmaxlagm0m0maxlag <- seq(-max.lag, max.lag)
}
if (m0 >= 1) {
seqmaxlagm0m0maxlag <- c(seq(-max.lag, -m0), seq(m0, max.lag))
}
seqm0maxlag <- seq(m0, max.lag)
cc <- results.tq[, 1]
t.tilde <- results.tq[, 2]
q.tilde.p <- results.tq[, 3]
q.tilde <- q.tilde.p[!is.na(q.tilde.p)]
t <- results.tq[, 4]
hb.p <- results.tq[, 5]
hb <- hb.p[!is.na(hb.p)]
z.cv <- qnorm(1 - alpha / 2)
vec1 <- rep(1, 2 * max.lag + 1)
s.cb <- z.cv / sqrt(n) * vec1
r.cb <- z.cv * cc / t.tilde
lu.s.cb <- cbind(-s.cb, s.cb)
lu.r.cb <- cbind(-r.cb, r.cb)
rownames(lu.r.cb) <- rownames(lu.s.cb)
colnames(lu.r.cb) <- colnames(lu.s.cb)
pv.t <- 2 * (1 - pnorm(abs(t)))
pv.ttilde <- 2 * (1 - pnorm(abs(t.tilde)))
pv.hb <- 1 - pchisq(hb, abs(seqmaxlagm0m0maxlag) - m0 + 1)
pv.qtilde <- 1 - pchisq(q.tilde, abs(seqmaxlagm0m0maxlag) - m0 + 1)
if (NROW(plot) == 1) {
plotC <- plot
plotQ <- plot
}
if (NROW(plot) == 2) {
plotC <- plot[1]
plotQ <- plot[2]
}
if (plotC == TRUE & max.lag >= 1) {
myplotcorr(max.lag, seqmaxlagmaxlag, cc, s.cb, r.cb, alpha, n, 2, my.names, scale.font)
}
if (plotQ == TRUE & max.lag >= m0 + 1) {
if (m0 == 0) {
df.chisq.cv <- abs(seqmaxlagmaxlag) + 1
}
if (m0 >= 1) {
df.chisq.cv <- c(rev(seqm0maxlag) - m0 + 1, rep(NA, 2 * m0 - 1), seqm0maxlag - m0 + 1)
}
chisq.cv <- qchisq(alpha, df.chisq.cv, lower.tail = FALSE)
Sys.sleep(2)
myplotstat(max.lag, seqmaxlagmaxlag, hb.p, q.tilde.p, alpha, chisq.cv, n, 2, my.names, " HB", expression(~tilde(Q)), 1, scale.font)
}
colnames(seqmaxlagmaxlag) <- NULL
rownames(seqmaxlagmaxlag) <- NULL
colnames(cc) <- NULL
rownames(cc) <- NULL
colnames(lu.s.cb) <- NULL
rownames(lu.s.cb) <- NULL
colnames(lu.r.cb) <- NULL
rownames(lu.r.cb) <- NULL
colnames(t) <- NULL
rownames(t) <- NULL
colnames(pv.t) <- NULL
rownames(pv.t) <- NULL
colnames(t.tilde) <- NULL
rownames(t.tilde) <- NULL
colnames(pv.ttilde) <- NULL
rownames(pv.ttilde) <- NULL
colnames(seqmaxlagm0m0maxlag) <- NULL
rownames(seqmaxlagm0m0maxlag) <- NULL
colnames(hb) <- NULL
rownames(hb) <- NULL
colnames(pv.hb) <- NULL
rownames(pv.hb) <- NULL
colnames(q.tilde) <- NULL
rownames(q.tilde) <- NULL
colnames(pv.qtilde) <- NULL
rownames(pv.qtilde) <- NULL
cc.test.out <- structure(list(lag = as.numeric(seqmaxlagmaxlag), cc = as.numeric(cc), scb = as.matrix(lu.s.cb), rcb = as.matrix(lu.r.cb), t = as.numeric(t), pvt = as.numeric(pv.t), ttilde = as.numeric(t.tilde), pvttilde = as.numeric(pv.ttilde), lagc = as.numeric(seqmaxlagm0m0maxlag), hb = as.numeric(hb), pvhb = as.numeric(pv.hb), qtilde = as.numeric(q.tilde), pvqtilde = as.numeric(pv.qtilde), alpha = alpha, varnames = my.names), class = "cc.test")
return(cc.test.out)
}
#' @export
print.cc.test <- function(x, ...) {
options(max.print = 1000000)
results.t <- cbind(x$t, x$pvt, x$ttilde, x$pvttilde)
results.t <- format(round(results.t, 3), nsmall = 3)
results.q <- cbind(x$hb, x$pvhb, x$qtilde, x$pvqtilde)
max.lag <- max(abs(x$lagc))
m0 <- min(abs(x$lagc))
if (m0 >= 1) {
results.q <- rbind(results.q[1:(max.lag - m0 + 1), ], matrix(NA, nrow = 2 * m0 - 1, ncol = 4), results.q[(max.lag - m0 + 2):(2 * (max.lag - m0 + 1)), ])
}
results.q <- format(round(results.q, 3), nsmall = 3)
results.cc <- x$cc
results.cc <- format(round(results.cc, 3), nsmall = 3)
results.cb.lu <- cbind(x$scb, x$rcb)
results.cb.lu <- format(round(results.cb.lu, 3), nsmall = 3)
results.scb <- paste("(", results.cb.lu[ , 1], ",", results.cb.lu[ , 2], ")", sep = "")
results.rcb <- paste("(", results.cb.lu[ , 3], ",", results.cb.lu[ , 4], ")", sep = "")
results.tq <- cbind(x$lag, results.cc, results.scb, results.rcb, x$lag, results.t, x$lag, results.q)
results.tq <- data.frame(results.tq, row.names = NULL)
names(results.tq) <- c("Lag", "CC", paste("Stand. CB(", 100 * (1 - x$alpha), "%)", sep = ""), paste("Robust CB(", 100 * (1 - x$alpha), "%)", sep = ""), "Lag", "t", "p-value", "t-tilde", "p-value", "Lag", "HB", "p-value", "Q-tilde", "p-value")
msg <- paste("Tests for zero cross-correlation of ", x$varnames[1], " and ", x$varnames[2] , "(-k)", sep = "")
table.tq <- kable(results.tq, align = "r", format = "pipe")
cat("\n")
cat(msg)
print(table.tq)
cat("\n")
}
cc.test.t.tmink <- function(x, y, max.lag, m0, alpha, lambda) {
n <- NROW(x)
x.tilde <- x - mean(x)
y.tilde <- y - mean(y)
matrix.e <- x.tilde * y.tilde
for (kk in 1:max.lag) {
matrix.e <- cbind(matrix.e, c(matrix(0, kk, 1), x.tilde[(kk + 1):n] * y.tilde[1:(n - kk)]))
}
matrix.esq <- matrix.e ^ 2
if (max.lag == 0) {
t.tilde <- sum(matrix.e) / sqrt(sum(matrix.esq))
}
if (max.lag >= 1) {
t.tilde <- colSums(matrix.e) / sqrt(colSums(matrix.esq))
}
big.rstar <- diag(max.lag + 1) / 2
if (max.lag >= 1) {
for (jj in 1:max.lag) {
for (kk in (jj + 1):(max.lag + 1)) {
max.jk <- max(jj, kk)
num.r <- t(matrix.e[, jj]) %*% matrix.e[, kk]
rjkstar.ind <- num.r / sqrt(t(matrix.esq[max.jk:n, jj]) %*% matrix.esq[max.jk:n, kk])
if (abs(rjkstar.ind) > lambda) {
big.rstar[jj, kk] <- num.r / sqrt(sum(matrix.esq[max.jk:n, jj]) * sum(matrix.esq[max.jk:n, kk]))
}
}
}
}
big.rstar <- big.rstar + t(big.rstar)
q.tilde <- matrix(NA, nrow = max.lag - m0 + 1, ncol = 1)
for (mm in (m0 + 1):(max.lag + 1)) {
q.tilde[mm - m0] <- t(t.tilde[(m0 + 1):mm]) %*% solve(big.rstar[(m0 + 1):mm, (m0 + 1):mm]) %*% t.tilde[(m0 + 1):mm]
}
cc <- ccf(x, y, max.lag, plot = FALSE)$acf[(max.lag + 1):(2 * max.lag + 1)]
t <- sqrt(n) * cc
seq0maxlag <- seq(0, max.lag)
sq.cc <- cc ^ 2
seq.overnmink <- (n - seq0maxlag) ^ (-1)
hb <- (n ^ 2) * cumsum(sq.cc[(m0 + 1):(max.lag + 1)] * seq.overnmink[(m0 + 1):(max.lag + 1)])
if (m0 == 0) {
return(cbind(cc, t.tilde, q.tilde, t, hb))
}
if (m0 >= 1) {
return(cbind(cc, t.tilde, c(rep(NA, m0), q.tilde), t, c(rep(NA, m0), hb)))
}
}
########################################################################################################################
### BIVARIATE STATISTICS: special case k=0
########################################################################################################################
#' Testing zero Pearson correlation
#'
#' The function rcorr.test computes the test statistics for examining the null hypothesis
#' of zero Pearson correlation for multivariate series in Dalla, Giraitis and Phillips (2022).
#'
#' @usage rcorr.test(x, plot = TRUE, var.names = NULL, scale.font = 1)
#' @param x A numeric matrix or a multivariate numeric time series object (ts, xts, zoo) or a data frame.
#' @param plot Logical. If TRUE the sample Pearson correlations and the p-values for significance are plotted. Default is TRUE.
#' @param var.names NULL or a character string specifying the variable names. If NULL and x has names, the names of x are used. If NULL and x has no names, the string c("x[1]","x[2]",...) is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The p-value of the robust \eqn{\widetilde{t}}{t-tilde} statistic is for testing the null hypothesis \eqn{H_0:\rho_{i,j}=0}{H[0]:\rho[i,j]=0},
#' where \eqn{\rho_{i,j}}{\rho[i,j]} denotes the correlation of \eqn{x_{i}}{x[i]} and \eqn{x_{j}}{x[j]}.
#' @return An object of class "rcorr.test", which is a list with the following components:
#' \item{pc}{The sample Pearson correlations.}
#' \item{pv}{The p-values for the \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{varnames}{The variable names used in the plot/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}. \cr
#' Giraitis, L., Li, Y. and Phillips, P. C. B. (2024). "Robust Inference on Correlation under General Heterogeneity". Journal of Econometrics, 244(1), 105691, \doi{doi:10.1016/j.jeconom.2024.105691}.
#' @examples
#' x <- matrix(rnorm(400), 100)
#' rcorr.test(x)
#' @importFrom stats cor pnorm is.ts
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom forcats fct_rev
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
rcorr.test <- function(x, plot = TRUE, var.names = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric matrix or numeric multivariate time series (ts, xts, zoo) object or numeric data frame')
}
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.names) & NROW(var.names) > 1 & !is.character(var.names)) stop('argument "var.names" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) == 1) stop('argument "x" must be multivariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (is.null(var.names) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop(paste('argument "x" must have', NCOL(x), 'names', sep = " "))
}
if (!is.null(var.names)) {
if (NROW(var.names) != NCOL(x)) stop(paste('argument "var.names" must contain', NCOL(x), 'names', sep = " "))
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
nv <- NCOL(x)
if (is.null(var.names)) {
if (!is.null(colnames(x))) {
my.names <- colnames(x)
label.xi <- FALSE
}
if (is.null(colnames(x))) {
my.names <- paste0("x", 1:nv)
label.xi <- TRUE
}
}
if (!is.null(var.names)) {
my.names <- var.names
label.xi <- FALSE
}
x <- as.matrix(x)
n <- NROW(x)
x.tilde <- x - t(colMeans(x) * matrix(1, nv, n))
pv <- matrix(0, nv, nv)
for (ii in 1:(nv - 1)) {
for (jj in (ii + 1):nv) {
t.tilde <- (t(x.tilde[, ii]) %*% x.tilde[, jj]) / sqrt(t(x.tilde[, ii] ^ 2) %*% (x.tilde[, jj] ^ 2))
pv[ii, jj] <- 2 * (1 - pnorm(abs(t.tilde)))
}
}
pv <- pv + t(pv)
diag(pv) <- NA
pc <- cor(x)
if (plot == TRUE) {
myplotcorrmat(pc, pv, my.names, scale.font, label.xi)
}
colnames(pc) <- NULL
rownames(pc) <- NULL
colnames(pv) <- NULL
rownames(pv) <- NULL
rcorr.test.out <- structure(list(pc = as.matrix(pc), pv = as.matrix(pv), varnames = my.names), class = "rcorr.test")
return(rcorr.test.out)
}
#' @export
print.rcorr.test <- function(x, ...) {
options(max.print = 1000000)
nv <- NCOL(x$pc)
df.pc <- as.data.frame(cbind(x$varnames, format(round(x$pc, 3), nsmall = 3)), stringsAsFactors = FALSE)
for (ii in 1:nv) {
df.pc[ii, (ii + 1)] <- 1
}
colnames(df.pc) <- c(" ", x$varnames)
df.pv <- as.data.frame(cbind(x$varnames, format(round(x$pv, 3), nsmall = 3)), stringsAsFactors = FALSE)
for (ii in 1:nv) {
df.pv[ii, ii + 1] <- " "
}
colnames(df.pv) <- c(" ", x$varnames)
cat("\n")
cat("Matrix of Pearson correlations\n")
cat("\n")
print(df.pc, row.names = FALSE)
cat("\n")
cat("Matrix of p-values\n")
cat("\n")
print(df.pv, row.names = FALSE)
cat("\n")
}
########################################################################################################################
### PLOT STATISTICS
########################################################################################################################
myplotcorr <- function(max.lag, seq.max.lag, ac.cc, s.cb, r.cb, alpha, n, uni.biv, my.names, scale.font) {
options(scipen = 999)
results <- cbind(seq.max.lag, ac.cc, s.cb, -s.cb, r.cb, -r.cb)
if (uni.biv == 1) {
colnames(results) <- c("lag", "AC", "scb", "scbl", "rcb", "rcbl")
}
if (uni.biv == 2) {
colnames(results) <- c("lag", "CC", "scb", "scbl", "rcb", "rcbl")
}
df.results <- data.frame(results)
my.ticks.x <- seq(seq.max.lag[1], max.lag, ceiling(0.1 * max.lag))
max.y.tick <- max(breaks_pretty()(c(0, results[, 2:6])))
min.y.tick <- min(breaks_pretty()(c(0, results[, 2:6])))
my.ticks.y <- breaks_pretty()(c(0, results[, 2:6]))
label.scb <- paste(" Standard CB(", 100 * (1 - alpha), "%)", sep = "")
label.rcb <- paste(" Robust CB(", 100 * (1 - alpha), "%)", sep = "")
if (scale.font != 1) {
scale.fontaxis <- 0.9 * scale.font
} else {
scale.fontaxis <- scale.font
}
scale.key <- 15 + 5 * scale.font
g.corr <- ggplot(df.results, aes(x = .data$lag)) +
theme_classic() +
theme(plot.title = element_text(size = 13 * scale.font, face = "bold"), legend.text = element_text(size = 12 * scale.font), axis.text = element_text(size = 10 * scale.fontaxis), axis.title = element_text(size = 10 * scale.font)) +
theme(legend.position = "top", legend.key.size = unit(scale.key, "pt"), legend.text = element_text(margin = margin(r = 10, unit = "pt"))) +
theme(axis.title.x = element_text(margin = margin(t = 10))) +
theme(axis.text.x = element_text(margin = margin(t = 4))) +
guides(fill = guide_legend(order = 1)) +
theme(axis.title.y = element_blank()) +
labs(x = "Lag") +
scale_x_continuous(breaks = my.ticks.x) +
scale_y_continuous(breaks = my.ticks.y, limits = c(min.y.tick, max.y.tick), expand = c(0, 0)) +
geom_blank(aes(y = max.y.tick)) +
geom_blank(aes(y = min.y.tick)) +
geom_hline(yintercept = 0)
if (uni.biv == 1) {
g.corr <- g.corr + geom_col(aes(y = .data$AC, fill = " AC"), width = 0.2) + scale_fill_manual("", breaks = " AC", values = "#5A8DDC")
}
if (uni.biv == 2) {
g.corr <- g.corr + geom_col(aes(y = .data$CC, fill = " CC"), width = 0.2) + scale_fill_manual("", breaks = " CC", values = "#5A8DDC")
}
g.corr <- g.corr + geom_line(aes(y = .data$scb, colour = "scb"), linetype = "dashed", linewidth = 1) +
geom_line(aes(y = .data$scbl), colour = "gray50", linetype = "dashed", linewidth = 1) +
geom_line(aes(y = .data$rcb, colour = "rcb"), linetype = "dashed", linewidth = 1) +
geom_line(aes(y = .data$rcbl), colour = "#B32026", linetype = "dashed", linewidth = 1) +
scale_colour_manual("", breaks = c("scb", "rcb"), values = c("gray50", "#B32026"), labels = c(label.scb, label.rcb))
if (uni.biv == 1) {
g.corr <- g.corr + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Autocorrelation of ", my.names, sep = ""))
}
if (uni.biv == 2) {
g.corr <- g.corr + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cross-correlation of ", my.names[1], " and ", my.names[2], "(-k)", sep = ""))
}
show(g.corr)
}
myplotstat <- function(max.lag, seq.max.lag, stat1.val, stat2.val, alpha, stat.cv, n, uni.biv, my.names, stat1.name, stat2.name, stat.id, scale.font) {
options(scipen = 999)
results <- cbind(seq.max.lag, stat1.val, stat2.val, stat.cv)
colnames(results) <- c("lag", "stat1", "stat2", "cv")
df.results <- data.frame(results)
my.ticks <- seq(seq.max.lag[1], max.lag, ceiling(0.1 * max.lag))
max.y.tick <- max(breaks_pretty()(c(0, results[, 2:4])))
min.y.tick <- min(breaks_pretty()(c(0, results[, 2:4])))
my.ticks.y <- breaks_pretty()(c(0, results[, 2:4]))
label.stat1 <- stat1.name
label.stat2 <- stat2.name
label.cv <- paste(" cv(", 100 * alpha, "%)", sep = "")
if (scale.font != 1) {
scale.fontaxis <- 0.9 * scale.font
} else {
scale.fontaxis <- scale.font
}
scale.key <- 15 + 5 * scale.font
g.stat <- ggplot(df.results, aes(x = .data$lag)) +
theme_classic() +
theme(plot.title = element_text(size = 13 * scale.font, face = "bold"), legend.text = element_text(size = 12 * scale.font), axis.text = element_text(size = 10 * scale.fontaxis), axis.title = element_text(size = 10 * scale.font)) +
theme(legend.position = "top", legend.key.size = unit(scale.key, "pt"), legend.text = element_text(margin = margin(r = 10, unit = "pt"))) +
theme(axis.title.x = element_text(margin = margin(t = 10))) +
theme(axis.text.x = element_text(margin = margin(t = 4))) +
guides(fill = guide_legend(order = 1)) +
theme(axis.title.y = element_blank()) +
labs(x = "Lag") +
scale_x_continuous(breaks = my.ticks) +
scale_y_continuous(breaks = my.ticks.y, limits = c(min.y.tick, max.y.tick), expand = c(0, 0)) +
geom_blank(aes(y = max.y.tick)) +
geom_blank(aes(y = min.y.tick)) +
geom_line(aes(y = .data$stat1, colour = "stat1", linetype = "stat1"), linewidth = 1, na.rm = TRUE) +
geom_line(aes(y = .data$stat2, colour = "stat2", linetype = "stat2"), linewidth = 1, na.rm = TRUE) +
geom_line(aes(y = .data$cv, colour = "cv", linetype = "cv"), linewidth = 1, na.rm = TRUE)
if (stat.id == 1) {
g.stat <- g.stat + scale_colour_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "gray50", "stat2" = "#B32026", "cv" = "black"), labels = c(label.stat1, label.stat2, label.cv)) +
scale_linetype_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "solid", "stat2" = "solid", "cv" = "dashed"), labels = c(label.stat1, label.stat2, label.cv))
}
if (stat.id == 2 | stat.id == 3) {
g.stat <- g.stat + scale_colour_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "#B32026", "stat2" = "gray50", "cv" = "black"), labels = c(label.stat1, label.stat2, label.cv)) +
scale_linetype_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "solid", "stat2" = "solid", "cv" = "dashed"), labels = c(label.stat1, label.stat2, label.cv))
}
if (uni.biv == 1 && stat.id == 1) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cumulative tests for zero autocorrelation of ", my.names, sep = ""))
}
if (uni.biv == 1 && stat.id == 2) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Tests for iid property of ", my.names, sep = ""))
}
if (uni.biv == 1 && stat.id == 3) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cumulative tests for iid property of ", my.names, sep = ""))
}
if (uni.biv == 2 && stat.id == 1) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cumulative tests for zero cross-correlation of ", my.names[1], " and ", my.names[2], "(-k)", sep = ""))
}
show(g.stat)
}
myplotcorrmat <- function(cc, pv, my.names, scale.font, label.xi) {
nv <- NCOL(cc)
colnames(pv) <- my.names
rownames(pv) <- my.names
melted_pv <- melt(pv)
colnames(melted_pv) <- c("Var1", "Var2", "pv")
colnames(cc) <- my.names
rownames(cc) <- my.names
melted_cc <- melt(cc)
colnames(melted_cc) <- c("Var1", "Var2", "cc")
ccpv <- cbind(melted_cc, melted_pv[, 3])
colnames(ccpv) <- c("Var1", "Var2", "cc", "pv")
labelcc <- format(round(ccpv[, 3], 3), nsmall = 3)
labelcc[is.na(ccpv[, 4])] <- 1
labelpv <- paste("(", format(round(ccpv[, 4], 3), nsmall = 3), ")", sep = "")
labelpv[is.na(ccpv[, 4])] <- ""
mylabel <- paste(labelcc, "\n", labelpv, sep = "")
if (label.xi == TRUE) {
labelticks <- matrix(NA, nv, 1)
for (ii in 1:nv) {
labelticks[ii] <- as.expression(bquote(x[.(ii)]))
}
} else {
labelticks <- my.names
}
ccpv[is.na(ccpv)] <- 1.3
mycolors <- c("#B32026", "#C85C61", "#DC989B", "#F1D5D6", "#FFFFFF", "#BEBEBE")
colcode <- cut(ccpv[, 4], breaks = c(0, 0.001, 0.01, 0.05, 0.1, 1), include.lowest = TRUE)
brk <- levels(colcode)
colcode2 <- cut(ccpv[, 4], breaks = c(0, 0.001, 0.01, 0.05, 0.1, 1, 1.5), include.lowest = TRUE)
mylabelspv <- c("[0,0.001]", "(0.001,0.01]", "(0.01,0.05]", "(0.05,0.1]", "(0.1,1]")
scale.key <- 5 + 10 * scale.font
g <- ggplot(ccpv, aes(.data$Var1, rev(.data$Var2))) +
geom_tile(aes(fill = colcode2), colour = "black", show.legend = TRUE) +
scale_fill_manual(name = "p-value", breaks = brk, values = mycolors, labels = mylabelspv, drop = FALSE) +
scale_x_discrete(position = "top", labels = labelticks) +
scale_y_discrete(labels = rev(labelticks)) +
theme(plot.title = element_text(size = 13 * scale.font, face = "bold"), legend.title = element_text(size = 12 * scale.font), legend.text = element_text(size = 11 * scale.font), axis.text = element_text(size = 11 * scale.font)) +
theme(legend.key.size = unit(scale.key, "pt")) +
theme(plot.title = element_text(hjust = 0.5)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(panel.background = element_blank(), axis.ticks = element_blank()) +
theme(axis.ticks.length = unit(0, "pt")) +
theme(axis.text.x = element_text(colour = "black"), axis.text.y = element_text(colour = "black")) +
theme(plot.title = element_text(margin = margin(b = 20))) +
geom_text(aes(.data$Var1, rev(.data$Var2), label = mylabel), color = "black", size = 3.5 * scale.font) +
ggtitle("Pearson correlations and p-values for significance")
show(g)
}
Try the testcorr package in your browser
Any scripts or data that you put into this service are public.
testcorr documentation built on April 20, 2026, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions iid.test print.ac.test ac.test
Documented in ac.test iid.test
########################################################################################################################
### UNIVARIATE STATISTICS
########################################################################################################################
#' Testing zero autocorrelation
#'
#' The function ac.test computes the test statistics for examining the null hypothesis
#' of zero autocorrelation for univariate time series given in Dalla, Giraitis and Phillips (2022).
#'
#' @usage
#' ac.test(x, max.lag, m0 = 1, alpha = 0.05, lambda = 2.576,
#' plot = TRUE, var.name = NULL, scale.font = 1)
#' @param x A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param max.lag Maximum lag at which to calculate the test statistics.
#' @param m0 Minimum lag at which to calculate the cumulative test statistics. Default is 1.
#' @param alpha Significance level for hypothesis testing used in the plots. Default is 0.05.
#' @param lambda Threshold in \eqn{\widetilde{Q}}{Q-tilde} test statistics. Default is 2.576.
#' @param plot Logical. If TRUE, 1) the sample autocorrelations with their confidence bands are plotted and 2) the cumulative test statistics with their critical values are plotted. Default is TRUE. Can be a logical vector for each of the plots 1)-2).
#' @param var.name NULL or a character string specifying the variable name. If NULL and x has name, the name of x is used. If NULL and x has no name, the string "x" is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The standard \eqn{t} and robust \eqn{\widetilde{t}}{t-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_k=0}{H[0]:\rho[k]=0} at lags \eqn{k=1,...,max.lag},
#' and the standard \eqn{LB} and robust \eqn{\widetilde{Q}}{Q-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_{m_0}=...=\rho_m=0}{H[0]:\rho[m0]=...=\rho[m]=0} at lags \eqn{m=m_0,...,max.lag}{m=m0,...,max.lag},
#' where \eqn{\rho_k}{\rho[k]} denotes the autocorrelation of \eqn{x_t}{x[t]} at lag \eqn{k}.
#' @return An object of class "ac.test", which is a list with the following components:
#' \item{lag}{The lags of the sample autocorrelations.}
#' \item{ac}{The sample autocorrelations.}
#' \item{scb}{The lower and upper limit of the confidence bands based on the standard test statistics.}
#' \item{rcb}{The lower and upper limit of the confidence bands based on the robust test statistics.}
#' \item{t}{The \eqn{t}{t} test statistics.}
#' \item{pvt}{The p-values for the \eqn{t}{t} test statistics.}
#' \item{ttilde}{The \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{pvttilde}{The p-values for the \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{lagc}{The lags of the cumulative test statistics.}
#' \item{lb}{The \eqn{LB} test statistics.}
#' \item{pvlb}{The p-values for the \eqn{LB} test statistics.}
#' \item{qtilde}{The \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{pvqtilde}{The p-values for the \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{alpha}{Significance level for hypothesis testing used in the plots.}
#' \item{varname}{The variable name used in the plots/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}. \cr
#' Giraitis, L., Li, Y. and Phillips, P. C. B. (2024). "Robust Inference on Correlation under General Heterogeneity". Journal of Econometrics, 244(1), 105691, \doi{doi:10.1016/j.jeconom.2024.105691}.
#' @examples
#' x <- rnorm(100)
#' ac.test(x, max.lag = 10)
#' @importFrom stats acf pchisq pnorm qchisq qnorm is.ts
#' @import ggplot2
#' @importFrom scales breaks_pretty
#' @importFrom knitr kable
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
ac.test <- function(x, max.lag, m0 = 1, alpha = 0.05, lambda = 2.576, plot = TRUE, var.name = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (!is.numeric(max.lag)) stop('argument "max.lag" must be numeric')
if (!is.numeric(m0)) stop('argument "m0" must be numeric')
if (!is.numeric(alpha)) stop('argument "alpha" must be numeric')
if (!is.numeric(lambda)) stop('argument "lambda" must be numeric')
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.name) & NROW(var.name) == 1 & !is.character(var.name)) stop('argument "var.name" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) != 1) stop('argument "x" must be univariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (max.lag != ceiling(max.lag) | max.lag < m0 | max.lag > (NROW(x) - 1)) stop('argument "max.lag" must be an integer value greater or equal to "m0" and less than the sample size')
if (m0 != ceiling(m0) | m0 < 1 | m0 > max.lag) stop('argument "m0" must be an integer value greater or equal to one and less or equal than "maxlag"')
if (alpha < 0 | alpha > 1) stop('argument "alpha" must be between 0 and 1')
if (lambda < 0) stop('argument "lambda" must be non-negative')
if (is.null(var.name) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop('argument "x" must have one name')
}
if (!is.null(var.name)) {
if (NROW(var.name) != NCOL(x)) stop('argument "var.name" must contain one name')
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
if (is.null(var.name)) {
if (!is.null(colnames(x))) {
my.name <- colnames(x)
}
if (is.null(colnames(x))) {
my.name <- "x"
}
}
if (!is.null(var.name)) {
my.name <- var.name
}
x <- as.numeric(x)
n <- NROW(x)
x.tilde <- x - mean(x)
seq1maxlag <- seq(1, max.lag)
seqm0maxlag <- seq(m0, max.lag)
matrix.e <- x.tilde[2:n] * x.tilde[1:(n - 1)]
if (max.lag >= 2) {
for (kk in 2:max.lag) {
matrix.e <- cbind(matrix.e, c(matrix(0, (kk - 1), 1), x.tilde[(kk + 1):n] * x.tilde[1:(n - kk)]))
}
}
matrix.esq <- matrix.e ^ 2
if (max.lag == 1) {
t.tilde <- sum(matrix.e) / sqrt(sum(matrix.esq))
}
if (max.lag >= 2) {
t.tilde <- colSums(matrix.e) / sqrt(colSums(matrix.esq))
}
big.rstar <- diag(max.lag) / 2
if (max.lag >= 2) {
for (jj in 1:(max.lag - 1)) {
for (kk in (jj + 1):max.lag) {
max.jk <- max(jj, kk)
num.r <- t(matrix.e[, jj]) %*% matrix.e[, kk]
rjkstar.ind <- num.r / sqrt(t(matrix.esq[max.jk:(n - 1), jj]) %*% matrix.esq[max.jk:(n - 1), kk])
if (abs(rjkstar.ind) > lambda) {
big.rstar[jj, kk] <- num.r / sqrt(sum(matrix.esq[max.jk:(n - 1), jj]) * sum(matrix.esq[max.jk:(n - 1), kk]))
}
}
}
}
big.rstar <- big.rstar + t(big.rstar)
q.tilde <- matrix(NA, nrow = max.lag - m0 + 1, ncol = 1)
for (mm in m0:max.lag) {
q.tilde[mm - m0 + 1] <- t(t.tilde[m0:mm]) %*% solve(big.rstar[m0:mm, m0:mm]) %*% t.tilde[m0:mm]
}
ac <- acf(x, max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
t <- sqrt(n) * ac
sq.ac <- ac ^ 2
seq.overnmink <- (n - seq1maxlag) ^ (-1)
lb <- n * (n + 2) * cumsum(sq.ac[m0:max.lag] * seq.overnmink[m0:max.lag])
z.cv <- qnorm(1 - alpha / 2)
s.cb <- z.cv / sqrt(n) * rep(1, max.lag)
r.cb <- z.cv * ac / t.tilde
lu.s.cb <- cbind(-s.cb, s.cb)
lu.r.cb <- cbind(-r.cb, r.cb)
rownames(lu.r.cb) <- rownames(lu.s.cb)
colnames(lu.r.cb) <- colnames(lu.s.cb)
pv.t <- 2 * (1 - pnorm(abs(t)))
pv.ttilde <- 2 * (1 - pnorm(abs(t.tilde)))
pv.lb <- 1 - pchisq(lb, seqm0maxlag - m0 + 1)
pv.qtilde <- 1 - pchisq(q.tilde, seqm0maxlag - m0 + 1)
if (NROW(plot) == 1) {
plotA <- plot
plotQ <- plot
}
if (NROW(plot) == 2) {
plotA <- plot[1]
plotQ <- plot[2]
}
if (plotA == TRUE & max.lag >= 2) {
myplotcorr(max.lag, seq1maxlag, ac, s.cb, r.cb, alpha, n, 1, my.name, scale.font)
}
if (plotQ == TRUE & max.lag >= m0 + 1) {
chisq.cv <- qchisq(alpha, seqm0maxlag - m0 + 1, lower.tail = FALSE)
Sys.sleep(2)
myplotstat(max.lag, seqm0maxlag, lb, q.tilde, alpha, chisq.cv, n, 1, my.name, " LB", expression(~tilde(Q)), 1, scale.font)
}
colnames(seq1maxlag) <- NULL
rownames(seq1maxlag) <- NULL
colnames(ac) <- NULL
rownames(ac) <- NULL
colnames(lu.s.cb) <- NULL
rownames(lu.s.cb) <- NULL
colnames(lu.r.cb) <- NULL
rownames(lu.r.cb) <- NULL
colnames(t) <- NULL
rownames(t) <- NULL
colnames(pv.t) <- NULL
rownames(pv.t) <- NULL
colnames(t.tilde) <- NULL
rownames(t.tilde) <- NULL
colnames(pv.ttilde) <- NULL
rownames(pv.ttilde) <- NULL
colnames(seqm0maxlag) <- NULL
rownames(seqm0maxlag) <- NULL
colnames(lb) <- NULL
rownames(lb) <- NULL
colnames(pv.lb) <- NULL
rownames(pv.lb) <- NULL
colnames(q.tilde) <- NULL
rownames(q.tilde) <- NULL
colnames(pv.qtilde) <- NULL
rownames(pv.qtilde) <- NULL
ac.test.out <- structure(list(lag = as.numeric(seq1maxlag), ac = as.numeric(ac), scb = as.matrix(lu.s.cb), rcb = as.matrix(lu.r.cb), t = as.numeric(t), pvt = as.numeric(pv.t), ttilde = as.numeric(t.tilde), pvttilde = as.numeric(pv.ttilde), lagc = as.numeric(seqm0maxlag), lb = as.numeric(lb), pvlb = as.numeric(pv.lb), qtilde = as.numeric(q.tilde), pvqtilde = as.numeric(pv.qtilde), alpha = alpha, varname = my.name), class = "ac.test")
return(ac.test.out)
}
#' @export
print.ac.test <- function(x, ...) {
options(max.print = 1000000)
results.t <- cbind(x$t, x$pvt, x$ttilde, x$pvttilde)
results.t <- format(round(results.t, 3), nsmall = 3)
m0 <- min(x$lagc)
if (m0 == 1) {
results.q <- cbind(x$lb, x$pvlb, x$qtilde, x$pvqtilde)
}
if (m0 >= 2) {
results.q <- rbind(matrix(NA, nrow = m0 - 1, ncol = 4), cbind(x$lb, x$pvlb, x$qtilde, x$pvqtilde))
}
results.q <- format(round(results.q, 3), nsmall = 3)
results.ac <- x$ac
results.ac <- format(round(results.ac, 3), nsmall = 3)
results.cb.lu <- cbind(x$scb, x$rcb)
results.cb.lu <- format(round(results.cb.lu, 3), nsmall = 3)
results.scb <- paste("(", results.cb.lu[ , 1], ",", results.cb.lu[ , 2], ")", sep = "")
results.rcb <- paste("(", results.cb.lu[ , 3], ",", results.cb.lu[ , 4], ")", sep = "")
results.tq <- cbind(x$lag, results.ac, results.scb, results.rcb, x$lag, results.t, x$lag, results.q)
results.tq <- data.frame(results.tq, row.names = NULL)
names(results.tq) <- c("Lag", "AC", paste("Stand. CB(", 100 * (1 - x$alpha), "%)", sep = ""), paste("Robust CB(", 100 * (1 - x$alpha), "%)", sep = ""), "Lag", "t", "p-value", "t-tilde", "p-value", "Lag", "LB", "p-value", "Q-tilde", "p-value")
msg <- paste("Tests for zero autocorrelation of ", x$varname, sep = "")
table.tq <- kable(results.tq, align = "r", format = "pipe")
cat("\n")
cat(msg)
print(table.tq)
cat("\n")
}
########################################################################################################################
### UNIVARIATE STATISTICS: iid
########################################################################################################################
#' Testing iid property
#'
#' The function iid.test computes the test statistics for examining the null hypothesis
#' of i.i.d. property for univariate series given in Dalla, Giraitis and Phillips (2022).
#'
#' @usage
#' iid.test(x, max.lag, m0 = 1, alpha = 0.05,
#' plot = TRUE, var.name = NULL, scale.font = 1)
#' @param x A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param max.lag Maximum lag at which to calculate the test statistics.
#' @param m0 Minimum lag at which to calculate the cumulative test statistics. Default is 1.
#' @param alpha Significance level for hypothesis testing used in the plots. Default is 0.05.
#' @param plot Logical. If TRUE, 1) the test statistics (J) and their critical values are plotted and 2) the cumulative test statistics (C) with their critical values are plotted. Default is TRUE. Can be a logical vector for each of the plots 1)-2).
#' @param var.name NULL or a character string specifying the variable name. If NULL and x has name, the name of x is used. If NULL and x has no name, the string "x" is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The \eqn{J_{x,|x|}}{J[x,|x|]} and \eqn{J_{x,x^2}}{J[x,x^2]} statistics are for testing the null hypothesis of i.i.d. at lag \eqn{k}, \eqn{k=1,...,max.lag},
#' and the \eqn{C_{x,|x|}}{C[x,|x|]} and \eqn{C_{x,x^2}}{C[x,x^2]} statistics are for testing the null hypothesis of i.i.d. at lags \eqn{m_0,...,m}, \eqn{m=m_0,...,max.lag}.
#' @return An object of class "iid.test", which is a list with the following components:
#' \item{lag}{The lags of the test statistics.}
#' \item{jab}{The \eqn{J_{x,|x|}}{J[x,|x|]} test statistics.}
#' \item{pvjab}{The p-values for the \eqn{J_{x,|x|}}{J[x,|x|]} test statistics.}
#' \item{jsq}{The \eqn{J_{x,x^2}}{J[x,x^2]} test statistics.}
#' \item{pvjsq}{The p-values for the \eqn{J_{x,x^2}}{J[x,x^2]} test statistics.}
#' \item{lagc}{The lags of the cumulative test statistics.}
#' \item{cab}{The \eqn{C_{x,|x|}}{C[x,|x|]} test statistics.}
#' \item{pvcab}{The p-values for the \eqn{C_{x,|x|}}{C[x,|x|]} test statistics.}
#' \item{csq}{The \eqn{C_{x,x^2}}{C[x,x^2]} test statistics.}
#' \item{pvcsq}{The p-values for the \eqn{C_{x,x^2}}{C[x,x^2]} test statistics.}
#' \item{alpha}{Significance level for hypothesis testing used in the plots.}
#' \item{varname}{The variable name used in the plots/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}.
#' @examples
#' x <- rnorm(100)
#' iid.test(x, max.lag = 10)
#' @importFrom stats acf pchisq pnorm qchisq qnorm is.ts
#' @import ggplot2
#' @importFrom scales breaks_pretty
#' @importFrom knitr kable
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
iid.test <- function(x, max.lag, m0 = 1, alpha = 0.05, plot = TRUE, var.name = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (!is.numeric(max.lag)) stop('argument "max.lag" must be numeric')
if (!is.numeric(m0)) stop('argument "m0" must be numeric')
if (!is.numeric(alpha)) stop('argument "alpha" must be numeric')
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.name) & NROW(var.name) == 1 & !is.character(var.name)) stop('argument "var.name" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) != 1) stop('argument "x" must be univariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (max.lag != ceiling(max.lag) | max.lag < m0 | max.lag > (NROW(x) - 1)) stop('argument "max.lag" must be an integer value greater or equal to "m0" and less than the sample size')
if (m0 != ceiling(m0) | m0 < 1 | m0 > max.lag) stop('argument "m0" must be an integer value greater or equal to one and less or equal than "maxlag"')
if (alpha < 0 | alpha > 1) stop('argument "alpha" must be between 0 and 1')
if (is.null(var.name) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop('argument "x" must have one name')
}
if (!is.null(var.name)) {
if (NROW(var.name) != NCOL(x)) stop('argument "var.name" must contain one name')
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
if (is.null(var.name)) {
if (!is.null(colnames(x))) {
my.name <- colnames(x)
}
if (is.null(colnames(x))) {
my.name <- "x"
}
}
if (!is.null(var.name)) {
my.name <- var.name
}
x <- as.numeric(x)
n <- NROW(x)
x.tilde <- x - mean(x)
seq1maxlag <- seq(1, max.lag)
seqm0maxlag <- seq(m0, max.lag)
ac <- acf(x, max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
ac.abs <- acf(abs(x.tilde), max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
ac.sq <- acf(x.tilde ^ 2, max.lag, plot = FALSE)$acf[2:(max.lag + 1)]
j.abs <- (n ^ 2) / (n - seq1maxlag) * (ac ^ 2 + ac.abs ^ 2)
j.sq <- (n ^ 2) / (n - seq1maxlag) * (ac ^ 2 + ac.sq ^ 2)
c.abs <- cumsum(j.abs[m0:max.lag])
c.sq <- cumsum(j.sq[m0:max.lag])
pv.jabs <- 1 - pchisq(j.abs, 2)
pv.jsq <- 1 - pchisq(j.sq, 2)
pv.cabs <- 1 - pchisq(c.abs, 2 * (seqm0maxlag - m0 + 1))
pv.csq <- 1 - pchisq(c.sq, 2 * (seqm0maxlag - m0 + 1))
if (NROW(plot) == 1) {
plotJ <- plot
plotC <- plot
}
if (NROW(plot) == 2) {
plotJ <- plot[1]
plotC <- plot[2]
}
if (plotJ == TRUE & max.lag >= 2) {
chisq2.cv <- rep(qchisq(alpha, 2, lower.tail = FALSE), max.lag)
fvar.name <- substr(my.name, 1, 1)
sub.abs <- bquote(.(fvar.name) * ",|" * .(fvar.name) * "|")
sub.sq <- bquote(.(fvar.name) * "," * .(fvar.name)^2)
name.j.abs <- bquote(~J[.(sub.abs)])
name.j.sq <- bquote(~J[.(sub.sq)])
myplotstat(max.lag, seq1maxlag, j.abs, j.sq, alpha, chisq2.cv, n, 1, my.name, name.j.abs, name.j.sq, 2, scale.font)
}
if (plotC == TRUE & max.lag >= m0 + 1) {
chisq.cv <- qchisq(alpha, 2 * (seqm0maxlag - m0 + 1), lower.tail = FALSE)
fvar.name <- substr(my.name, 1, 1)
sub.abs <- bquote(.(fvar.name) * ",|" * .(fvar.name) * "|")
sub.sq <- bquote(.(fvar.name) * "," * .(fvar.name)^2)
name.c.abs <- bquote(~C[.(sub.abs)])
name.c.sq <- bquote(~C[.(sub.sq)])
Sys.sleep(2)
myplotstat(max.lag, seqm0maxlag, c.abs, c.sq, alpha, chisq.cv, n, 1, my.name, name.c.abs, name.c.sq, 3, scale.font)
}
colnames(seq1maxlag) <- NULL
rownames(seq1maxlag) <- NULL
colnames(j.abs) <- NULL
rownames(j.abs) <- NULL
colnames(pv.jabs) <- NULL
rownames(pv.jabs) <- NULL
colnames(j.sq) <- NULL
rownames(j.sq) <- NULL
colnames(pv.jsq) <- NULL
rownames(pv.jsq) <- NULL
colnames(seqm0maxlag) <- NULL
rownames(seqm0maxlag) <- NULL
colnames(c.abs) <- NULL
rownames(c.abs) <- NULL
colnames(pv.cabs) <- NULL
rownames(pv.cabs) <- NULL
colnames(c.sq) <- NULL
rownames(c.sq) <- NULL
colnames(pv.csq) <- NULL
rownames(pv.csq) <- NULL
iid.test.out <- structure(list(lag = as.numeric(seq1maxlag), jabs = as.numeric(j.abs), pvjabs = as.numeric(pv.jabs), jsq = as.numeric(j.sq), pvjsq = as.numeric(pv.jsq), lagc = as.numeric(seqm0maxlag), cabs = as.numeric(c.abs), pvcabs = as.numeric(pv.cabs), csq = as.numeric(c.sq), pvcsq = as.numeric(pv.csq), alpha = alpha, varname = my.name), class = "iid.test")
return(iid.test.out)
}
#' @export
print.iid.test <- function(x, ...) {
options(max.print = 1000000)
results.j <- cbind(x$jabs, x$pvjabs, x$jsq, x$pvjsq)
results.j <- format(round(results.j, 3), nsmall = 3)
m0 <- min(x$lagc)
if (m0 == 1) {
results.c <- cbind(x$cabs, x$pvcabs, x$csq, x$pvcsq)
}
if (m0 >= 2) {
results.c <- rbind(matrix(NA, nrow = m0 - 1, ncol = 4), cbind(x$cabs, x$pvcabs, x$csq, x$pvcsq))
}
results.c <- format(round(results.c, 3), nsmall = 3)
results.jc <- cbind(x$lag, results.j, x$lag, results.c)
results.jc <- data.frame(results.jc, row.names = NULL)
fvar.name <- substr(x$varname, 1, 1)
name.j.abs <- paste("J[", fvar.name, ",", "\u01C0", fvar.name, "\u01C0]", sep = "")
name.j.sq <- paste("J[", fvar.name, ",", fvar.name, "\u00B2]", sep = "")
name.c.abs <- paste("C[", fvar.name, ",", "\u01C0", fvar.name, "\u01C0]", sep = "")
name.c.sq <- paste("C[", fvar.name, ",", fvar.name, "\u00B2]", sep = "")
names(results.jc) <- c("Lag", name.j.abs, "p-value", name.j.sq, "p-value", "Lag", name.c.abs, "p-value", name.c.sq, "p-value")
msg <- paste("Tests for i.i.d. property of ", x$varname, sep = "")
table.jc <- kable(results.jc, align = "r", format = "pipe")
cat("\n")
cat(msg)
print(table.jc)
cat("\n")
}
########################################################################################################################
### BIVARIATE STATISTICS
########################################################################################################################
#' Testing zero cross-correlation
#'
#' The function cc.test computes the test statistics for examining the null hypothesis
#' of zero cross-correlation for bivariate time series given in Dalla, Giraitis and Phillips (2022).
#'
#' @usage
#' cc.test(x, y, max.lag, m0 = 0, alpha = 0.05, lambda = 2.576,
#' plot = TRUE, var.names = NULL, scale.font = 1)
#' @param x A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param y A numeric vector or a univariate numeric time series (ts, xts, zoo) object or a data frame variable.
#' @param max.lag Maximum lag at which to calculate the test statistics.
#' @param m0 Minimum lag at which to calculate the cumulative test statistics. Default is 0.
#' @param alpha Significance level for hypothesis testing used in the plots. Default is 0.05.
#' @param lambda Threshold in \eqn{\widetilde{Q}}{Q-tilde} test statistics. Default is 2.576.
#' @param plot Logical. If TRUE, 1) the sample cross-correlations with their confidence bands are plotted and 2) the cumulative test statistics with their critical values are plotted. Default is TRUE. Can be a logical vector for each of the plots 1)-2).
#' @param var.names NULL or a character string specifying the variable names. If NULL and x,y have names, the names of x,y are used. If NULL and x,y have no names, the string c("x","y") is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The standard \eqn{t} and robust \eqn{\widetilde{t}}{t-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_k=0}{H[0]:\rho[k]=0} at lags \eqn{k=-max.lag,...,-1,0,1,max.lag},
#' and the standard \eqn{HB} and robust \eqn{\widetilde{Q}}{Q-tilde} statistics are for testing the null hypothesis \eqn{H_0:\rho_{m_0}=...=\rho_m=0}{H[0]:\rho[m0]=...=\rho[m]=0} at lags \eqn{m=-max.lag,...,-1,0,1,max.lag},
#' where \eqn{\rho_k}{\rho[k]} denotes the cross-correlation of \eqn{x_t}{x[t]} and \eqn{y_{t-k}}{y[t-k]} at lag \eqn{k}.
#' @return An object of class "cc.test", which is a list with the following components:
#' \item{lag}{The lags of the sample cross-correlations.}
#' \item{cc}{The sample cross-correlations.}
#' \item{scb}{The lower and upper limit of the confidence bands based on the standard test statistics.}
#' \item{rcb}{The lower and upper limit of the confidence bands based on the robust test statistics.}
#' \item{t}{The \eqn{t}{t} test statistics.}
#' \item{pvt}{The p-values for the \eqn{t}{t} test statistics.}
#' \item{ttilde}{The \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{pvtttilde}{The p-values for the \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{lagc}{The lags of the cumulative test statistics.}
#' \item{hb}{The \eqn{HB} test statistics.}
#' \item{pvhb}{The p-values for the \eqn{HB} test statistics.}
#' \item{qtilde}{The \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{pvqtilde}{The p-values for the \eqn{\widetilde{Q}}{Q-tilde} test statistics.}
#' \item{alpha}{Significance level for hypothesis testing used in the plots.}
#' \item{varnames}{The variable names used in the plots/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}. \cr
#' Giraitis, L., Li, Y. and Phillips, P. C. B. (2024). "Robust Inference on Correlation under General Heterogeneity". Journal of Econometrics, 244(1), 105691, \doi{doi:10.1016/j.jeconom.2024.105691}.
#' @examples
#' x <- rnorm(100)
#' y <- rnorm(100)
#' cc.test(x, y, max.lag = 10)
#' @importFrom stats acf ccf pchisq pnorm qchisq qnorm is.ts
#' @import ggplot2
#' @importFrom scales breaks_pretty
#' @importFrom knitr kable
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
cc.test <- function(x, y, max.lag, m0 = 0, alpha = 0.05, lambda = 2.576, plot = TRUE, var.names = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.null(y)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (is.ts(y) | is.xts(y) | is.zoo(y) | is.numeric(y) | is.data.frame(y)) {}
else{
stop('argument "y" must be numeric vector or numeric time series (ts, xts, zoo) object or numeric data frame variable')
}
if (!is.numeric(max.lag)) stop('argument "max.lag" must be numeric')
if (!is.numeric(m0)) stop('argument "m0" must be numeric')
if (!is.numeric(alpha)) stop('argument "alpha" must be numeric')
if (!is.numeric(lambda)) stop('argument "lambda" must be numeric')
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.names) & NROW(var.names) == 2 & (!is.character(var.names[1]) | !is.character(var.names[2]))) stop('argument "var.names" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) != 1) stop('argument "x" must be univariate')
if (NCOL(y) != 1) stop('argument "y" must be univariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (any(is.na(y))) stop('argument "y" must not contain missing values')
if (NCOL(x) != NCOL(y)) stop('arguments "x" and "y" must have the same length')
if (max.lag != ceiling(max.lag) | max.lag < m0 | max.lag > (NROW(x) - 1)) stop('argument "max.lag" must be an integer value greater or equal to "m0" and less than the sample size')
if (m0 != ceiling(m0) | m0 < 0 | m0 > max.lag) stop('argument "m0" must be an integer value greater or equal to zero and less or equal than "maxlag"')
if (alpha < 0 | alpha > 1) stop('argument "alpha" must be between 0 and 1')
if (lambda < 0) stop('argument "lambda" must be non-negative')
if (is.null(var.names) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop('argument "x" must have one name')
}
if (is.null(var.names) & !is.null(colnames(y))) {
if (NROW(colnames(y)) != NCOL(x)) stop('argument "y" must have one name')
}
if (!is.null(var.names)) {
if (NROW(var.names) != 2) stop('argument "var.names" must contain two names')
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
if (is.null(var.names)) {
if (!is.null(colnames(x))) {
my.names <- colnames(x)
}
if (is.null(colnames(x))) {
my.names <- c("x","y")
}
}
if (!is.null(var.names)) {
my.names <- var.names
}
x <- as.numeric(x)
y <- as.numeric(y)
n <- NROW(x)
results.tq.xyk <- cc.test.t.tmink(x, y, max.lag, m0, alpha, lambda)
results.tq.yxk <- cc.test.t.tmink(y, x, max.lag, m0, alpha, lambda)
results.tq <- rbind(results.tq.yxk[(max.lag + 1):2, ], results.tq.xyk)
seqmaxlagmaxlag <- seq(-max.lag, max.lag)
if (m0 == 0) {
seqmaxlagm0m0maxlag <- seq(-max.lag, max.lag)
}
if (m0 >= 1) {
seqmaxlagm0m0maxlag <- c(seq(-max.lag, -m0), seq(m0, max.lag))
}
seqm0maxlag <- seq(m0, max.lag)
cc <- results.tq[, 1]
t.tilde <- results.tq[, 2]
q.tilde.p <- results.tq[, 3]
q.tilde <- q.tilde.p[!is.na(q.tilde.p)]
t <- results.tq[, 4]
hb.p <- results.tq[, 5]
hb <- hb.p[!is.na(hb.p)]
z.cv <- qnorm(1 - alpha / 2)
vec1 <- rep(1, 2 * max.lag + 1)
s.cb <- z.cv / sqrt(n) * vec1
r.cb <- z.cv * cc / t.tilde
lu.s.cb <- cbind(-s.cb, s.cb)
lu.r.cb <- cbind(-r.cb, r.cb)
rownames(lu.r.cb) <- rownames(lu.s.cb)
colnames(lu.r.cb) <- colnames(lu.s.cb)
pv.t <- 2 * (1 - pnorm(abs(t)))
pv.ttilde <- 2 * (1 - pnorm(abs(t.tilde)))
pv.hb <- 1 - pchisq(hb, abs(seqmaxlagm0m0maxlag) - m0 + 1)
pv.qtilde <- 1 - pchisq(q.tilde, abs(seqmaxlagm0m0maxlag) - m0 + 1)
if (NROW(plot) == 1) {
plotC <- plot
plotQ <- plot
}
if (NROW(plot) == 2) {
plotC <- plot[1]
plotQ <- plot[2]
}
if (plotC == TRUE & max.lag >= 1) {
myplotcorr(max.lag, seqmaxlagmaxlag, cc, s.cb, r.cb, alpha, n, 2, my.names, scale.font)
}
if (plotQ == TRUE & max.lag >= m0 + 1) {
if (m0 == 0) {
df.chisq.cv <- abs(seqmaxlagmaxlag) + 1
}
if (m0 >= 1) {
df.chisq.cv <- c(rev(seqm0maxlag) - m0 + 1, rep(NA, 2 * m0 - 1), seqm0maxlag - m0 + 1)
}
chisq.cv <- qchisq(alpha, df.chisq.cv, lower.tail = FALSE)
Sys.sleep(2)
myplotstat(max.lag, seqmaxlagmaxlag, hb.p, q.tilde.p, alpha, chisq.cv, n, 2, my.names, " HB", expression(~tilde(Q)), 1, scale.font)
}
colnames(seqmaxlagmaxlag) <- NULL
rownames(seqmaxlagmaxlag) <- NULL
colnames(cc) <- NULL
rownames(cc) <- NULL
colnames(lu.s.cb) <- NULL
rownames(lu.s.cb) <- NULL
colnames(lu.r.cb) <- NULL
rownames(lu.r.cb) <- NULL
colnames(t) <- NULL
rownames(t) <- NULL
colnames(pv.t) <- NULL
rownames(pv.t) <- NULL
colnames(t.tilde) <- NULL
rownames(t.tilde) <- NULL
colnames(pv.ttilde) <- NULL
rownames(pv.ttilde) <- NULL
colnames(seqmaxlagm0m0maxlag) <- NULL
rownames(seqmaxlagm0m0maxlag) <- NULL
colnames(hb) <- NULL
rownames(hb) <- NULL
colnames(pv.hb) <- NULL
rownames(pv.hb) <- NULL
colnames(q.tilde) <- NULL
rownames(q.tilde) <- NULL
colnames(pv.qtilde) <- NULL
rownames(pv.qtilde) <- NULL
cc.test.out <- structure(list(lag = as.numeric(seqmaxlagmaxlag), cc = as.numeric(cc), scb = as.matrix(lu.s.cb), rcb = as.matrix(lu.r.cb), t = as.numeric(t), pvt = as.numeric(pv.t), ttilde = as.numeric(t.tilde), pvttilde = as.numeric(pv.ttilde), lagc = as.numeric(seqmaxlagm0m0maxlag), hb = as.numeric(hb), pvhb = as.numeric(pv.hb), qtilde = as.numeric(q.tilde), pvqtilde = as.numeric(pv.qtilde), alpha = alpha, varnames = my.names), class = "cc.test")
return(cc.test.out)
}
#' @export
print.cc.test <- function(x, ...) {
options(max.print = 1000000)
results.t <- cbind(x$t, x$pvt, x$ttilde, x$pvttilde)
results.t <- format(round(results.t, 3), nsmall = 3)
results.q <- cbind(x$hb, x$pvhb, x$qtilde, x$pvqtilde)
max.lag <- max(abs(x$lagc))
m0 <- min(abs(x$lagc))
if (m0 >= 1) {
results.q <- rbind(results.q[1:(max.lag - m0 + 1), ], matrix(NA, nrow = 2 * m0 - 1, ncol = 4), results.q[(max.lag - m0 + 2):(2 * (max.lag - m0 + 1)), ])
}
results.q <- format(round(results.q, 3), nsmall = 3)
results.cc <- x$cc
results.cc <- format(round(results.cc, 3), nsmall = 3)
results.cb.lu <- cbind(x$scb, x$rcb)
results.cb.lu <- format(round(results.cb.lu, 3), nsmall = 3)
results.scb <- paste("(", results.cb.lu[ , 1], ",", results.cb.lu[ , 2], ")", sep = "")
results.rcb <- paste("(", results.cb.lu[ , 3], ",", results.cb.lu[ , 4], ")", sep = "")
results.tq <- cbind(x$lag, results.cc, results.scb, results.rcb, x$lag, results.t, x$lag, results.q)
results.tq <- data.frame(results.tq, row.names = NULL)
names(results.tq) <- c("Lag", "CC", paste("Stand. CB(", 100 * (1 - x$alpha), "%)", sep = ""), paste("Robust CB(", 100 * (1 - x$alpha), "%)", sep = ""), "Lag", "t", "p-value", "t-tilde", "p-value", "Lag", "HB", "p-value", "Q-tilde", "p-value")
msg <- paste("Tests for zero cross-correlation of ", x$varnames[1], " and ", x$varnames[2] , "(-k)", sep = "")
table.tq <- kable(results.tq, align = "r", format = "pipe")
cat("\n")
cat(msg)
print(table.tq)
cat("\n")
}
cc.test.t.tmink <- function(x, y, max.lag, m0, alpha, lambda) {
n <- NROW(x)
x.tilde <- x - mean(x)
y.tilde <- y - mean(y)
matrix.e <- x.tilde * y.tilde
for (kk in 1:max.lag) {
matrix.e <- cbind(matrix.e, c(matrix(0, kk, 1), x.tilde[(kk + 1):n] * y.tilde[1:(n - kk)]))
}
matrix.esq <- matrix.e ^ 2
if (max.lag == 0) {
t.tilde <- sum(matrix.e) / sqrt(sum(matrix.esq))
}
if (max.lag >= 1) {
t.tilde <- colSums(matrix.e) / sqrt(colSums(matrix.esq))
}
big.rstar <- diag(max.lag + 1) / 2
if (max.lag >= 1) {
for (jj in 1:max.lag) {
for (kk in (jj + 1):(max.lag + 1)) {
max.jk <- max(jj, kk)
num.r <- t(matrix.e[, jj]) %*% matrix.e[, kk]
rjkstar.ind <- num.r / sqrt(t(matrix.esq[max.jk:n, jj]) %*% matrix.esq[max.jk:n, kk])
if (abs(rjkstar.ind) > lambda) {
big.rstar[jj, kk] <- num.r / sqrt(sum(matrix.esq[max.jk:n, jj]) * sum(matrix.esq[max.jk:n, kk]))
}
}
}
}
big.rstar <- big.rstar + t(big.rstar)
q.tilde <- matrix(NA, nrow = max.lag - m0 + 1, ncol = 1)
for (mm in (m0 + 1):(max.lag + 1)) {
q.tilde[mm - m0] <- t(t.tilde[(m0 + 1):mm]) %*% solve(big.rstar[(m0 + 1):mm, (m0 + 1):mm]) %*% t.tilde[(m0 + 1):mm]
}
cc <- ccf(x, y, max.lag, plot = FALSE)$acf[(max.lag + 1):(2 * max.lag + 1)]
t <- sqrt(n) * cc
seq0maxlag <- seq(0, max.lag)
sq.cc <- cc ^ 2
seq.overnmink <- (n - seq0maxlag) ^ (-1)
hb <- (n ^ 2) * cumsum(sq.cc[(m0 + 1):(max.lag + 1)] * seq.overnmink[(m0 + 1):(max.lag + 1)])
if (m0 == 0) {
return(cbind(cc, t.tilde, q.tilde, t, hb))
}
if (m0 >= 1) {
return(cbind(cc, t.tilde, c(rep(NA, m0), q.tilde), t, c(rep(NA, m0), hb)))
}
}
########################################################################################################################
### BIVARIATE STATISTICS: special case k=0
########################################################################################################################
#' Testing zero Pearson correlation
#'
#' The function rcorr.test computes the test statistics for examining the null hypothesis
#' of zero Pearson correlation for multivariate series in Dalla, Giraitis and Phillips (2022).
#'
#' @usage rcorr.test(x, plot = TRUE, var.names = NULL, scale.font = 1)
#' @param x A numeric matrix or a multivariate numeric time series object (ts, xts, zoo) or a data frame.
#' @param plot Logical. If TRUE the sample Pearson correlations and the p-values for significance are plotted. Default is TRUE.
#' @param var.names NULL or a character string specifying the variable names. If NULL and x has names, the names of x are used. If NULL and x has no names, the string c("x[1]","x[2]",...) is used. Default is NULL.
#' @param scale.font A positive number indicating the scaling of the font size in the plots. Default is 1.
#' @details
#' The p-value of the robust \eqn{\widetilde{t}}{t-tilde} statistic is for testing the null hypothesis \eqn{H_0:\rho_{i,j}=0}{H[0]:\rho[i,j]=0},
#' where \eqn{\rho_{i,j}}{\rho[i,j]} denotes the correlation of \eqn{x_{i}}{x[i]} and \eqn{x_{j}}{x[j]}.
#' @return An object of class "rcorr.test", which is a list with the following components:
#' \item{pc}{The sample Pearson correlations.}
#' \item{pv}{The p-values for the \eqn{\widetilde{t}}{t-tilde} test statistics.}
#' \item{varnames}{The variable names used in the plot/table.}
#' @note
#' Missing values are not allowed.
#' @author
#' Violetta Dalla, Liudas Giraitis and Peter C. B. Phillips
#' @references
#' Dalla, V., Giraitis, L. and Phillips, P. C. B. (2022). "Robust Tests for White Noise and Cross-Correlation". Econometric Theory, 38(5), 913-941, \doi{doi:10.1017/S0266466620000341}. Cowles Foundation, Discussion Paper No. 2194RS, \url{https://elischolar.library.yale.edu/cowles-discussion-paper-series/57/}. \cr
#' Giraitis, L., Li, Y. and Phillips, P. C. B. (2024). "Robust Inference on Correlation under General Heterogeneity". Journal of Econometrics, 244(1), 105691, \doi{doi:10.1016/j.jeconom.2024.105691}.
#' @examples
#' x <- matrix(rnorm(400), 100)
#' rcorr.test(x)
#' @importFrom stats cor pnorm is.ts
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom forcats fct_rev
#' @importFrom methods show
#' @importFrom xts is.xts
#' @importFrom zoo is.zoo
#' @export
rcorr.test <- function(x, plot = TRUE, var.names = NULL, scale.font = 1) {
if (is.null(x)) stop()
if (is.ts(x) | is.xts(x) | is.zoo(x) | is.numeric(x) | is.data.frame(x)) {}
else{
stop('argument "x" must be numeric matrix or numeric multivariate time series (ts, xts, zoo) object or numeric data frame')
}
if (!is.logical(plot)) stop('argument "plot" must be logical')
if (!is.null(var.names) & NROW(var.names) > 1 & !is.character(var.names)) stop('argument "var.names" must be NULL or string')
if (!is.numeric(scale.font)) stop('argument "scale.font" must be numeric')
if (NCOL(x) == 1) stop('argument "x" must be multivariate')
if (any(is.na(x))) stop('argument "x" must not contain missing values')
if (is.null(var.names) & !is.null(colnames(x))) {
if (NROW(colnames(x)) != NCOL(x)) stop(paste('argument "x" must have', NCOL(x), 'names', sep = " "))
}
if (!is.null(var.names)) {
if (NROW(var.names) != NCOL(x)) stop(paste('argument "var.names" must contain', NCOL(x), 'names', sep = " "))
}
if (scale.font <= 0) stop('argument "scale.font" must be positive')
nv <- NCOL(x)
if (is.null(var.names)) {
if (!is.null(colnames(x))) {
my.names <- colnames(x)
label.xi <- FALSE
}
if (is.null(colnames(x))) {
my.names <- paste0("x", 1:nv)
label.xi <- TRUE
}
}
if (!is.null(var.names)) {
my.names <- var.names
label.xi <- FALSE
}
x <- as.matrix(x)
n <- NROW(x)
x.tilde <- x - t(colMeans(x) * matrix(1, nv, n))
pv <- matrix(0, nv, nv)
for (ii in 1:(nv - 1)) {
for (jj in (ii + 1):nv) {
t.tilde <- (t(x.tilde[, ii]) %*% x.tilde[, jj]) / sqrt(t(x.tilde[, ii] ^ 2) %*% (x.tilde[, jj] ^ 2))
pv[ii, jj] <- 2 * (1 - pnorm(abs(t.tilde)))
}
}
pv <- pv + t(pv)
diag(pv) <- NA
pc <- cor(x)
if (plot == TRUE) {
myplotcorrmat(pc, pv, my.names, scale.font, label.xi)
}
colnames(pc) <- NULL
rownames(pc) <- NULL
colnames(pv) <- NULL
rownames(pv) <- NULL
rcorr.test.out <- structure(list(pc = as.matrix(pc), pv = as.matrix(pv), varnames = my.names), class = "rcorr.test")
return(rcorr.test.out)
}
#' @export
print.rcorr.test <- function(x, ...) {
options(max.print = 1000000)
nv <- NCOL(x$pc)
df.pc <- as.data.frame(cbind(x$varnames, format(round(x$pc, 3), nsmall = 3)), stringsAsFactors = FALSE)
for (ii in 1:nv) {
df.pc[ii, (ii + 1)] <- 1
}
colnames(df.pc) <- c(" ", x$varnames)
df.pv <- as.data.frame(cbind(x$varnames, format(round(x$pv, 3), nsmall = 3)), stringsAsFactors = FALSE)
for (ii in 1:nv) {
df.pv[ii, ii + 1] <- " "
}
colnames(df.pv) <- c(" ", x$varnames)
cat("\n")
cat("Matrix of Pearson correlations\n")
cat("\n")
print(df.pc, row.names = FALSE)
cat("\n")
cat("Matrix of p-values\n")
cat("\n")
print(df.pv, row.names = FALSE)
cat("\n")
}
########################################################################################################################
### PLOT STATISTICS
########################################################################################################################
myplotcorr <- function(max.lag, seq.max.lag, ac.cc, s.cb, r.cb, alpha, n, uni.biv, my.names, scale.font) {
options(scipen = 999)
results <- cbind(seq.max.lag, ac.cc, s.cb, -s.cb, r.cb, -r.cb)
if (uni.biv == 1) {
colnames(results) <- c("lag", "AC", "scb", "scbl", "rcb", "rcbl")
}
if (uni.biv == 2) {
colnames(results) <- c("lag", "CC", "scb", "scbl", "rcb", "rcbl")
}
df.results <- data.frame(results)
my.ticks.x <- seq(seq.max.lag[1], max.lag, ceiling(0.1 * max.lag))
max.y.tick <- max(breaks_pretty()(c(0, results[, 2:6])))
min.y.tick <- min(breaks_pretty()(c(0, results[, 2:6])))
my.ticks.y <- breaks_pretty()(c(0, results[, 2:6]))
label.scb <- paste(" Standard CB(", 100 * (1 - alpha), "%)", sep = "")
label.rcb <- paste(" Robust CB(", 100 * (1 - alpha), "%)", sep = "")
if (scale.font != 1) {
scale.fontaxis <- 0.9 * scale.font
} else {
scale.fontaxis <- scale.font
}
scale.key <- 15 + 5 * scale.font
g.corr <- ggplot(df.results, aes(x = .data$lag)) +
theme_classic() +
theme(plot.title = element_text(size = 13 * scale.font, face = "bold"), legend.text = element_text(size = 12 * scale.font), axis.text = element_text(size = 10 * scale.fontaxis), axis.title = element_text(size = 10 * scale.font)) +
theme(legend.position = "top", legend.key.size = unit(scale.key, "pt"), legend.text = element_text(margin = margin(r = 10, unit = "pt"))) +
theme(axis.title.x = element_text(margin = margin(t = 10))) +
theme(axis.text.x = element_text(margin = margin(t = 4))) +
guides(fill = guide_legend(order = 1)) +
theme(axis.title.y = element_blank()) +
labs(x = "Lag") +
scale_x_continuous(breaks = my.ticks.x) +
scale_y_continuous(breaks = my.ticks.y, limits = c(min.y.tick, max.y.tick), expand = c(0, 0)) +
geom_blank(aes(y = max.y.tick)) +
geom_blank(aes(y = min.y.tick)) +
geom_hline(yintercept = 0)
if (uni.biv == 1) {
g.corr <- g.corr + geom_col(aes(y = .data$AC, fill = " AC"), width = 0.2) + scale_fill_manual("", breaks = " AC", values = "#5A8DDC")
}
if (uni.biv == 2) {
g.corr <- g.corr + geom_col(aes(y = .data$CC, fill = " CC"), width = 0.2) + scale_fill_manual("", breaks = " CC", values = "#5A8DDC")
}
g.corr <- g.corr + geom_line(aes(y = .data$scb, colour = "scb"), linetype = "dashed", linewidth = 1) +
geom_line(aes(y = .data$scbl), colour = "gray50", linetype = "dashed", linewidth = 1) +
geom_line(aes(y = .data$rcb, colour = "rcb"), linetype = "dashed", linewidth = 1) +
geom_line(aes(y = .data$rcbl), colour = "#B32026", linetype = "dashed", linewidth = 1) +
scale_colour_manual("", breaks = c("scb", "rcb"), values = c("gray50", "#B32026"), labels = c(label.scb, label.rcb))
if (uni.biv == 1) {
g.corr <- g.corr + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Autocorrelation of ", my.names, sep = ""))
}
if (uni.biv == 2) {
g.corr <- g.corr + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cross-correlation of ", my.names[1], " and ", my.names[2], "(-k)", sep = ""))
}
show(g.corr)
}
myplotstat <- function(max.lag, seq.max.lag, stat1.val, stat2.val, alpha, stat.cv, n, uni.biv, my.names, stat1.name, stat2.name, stat.id, scale.font) {
options(scipen = 999)
results <- cbind(seq.max.lag, stat1.val, stat2.val, stat.cv)
colnames(results) <- c("lag", "stat1", "stat2", "cv")
df.results <- data.frame(results)
my.ticks <- seq(seq.max.lag[1], max.lag, ceiling(0.1 * max.lag))
max.y.tick <- max(breaks_pretty()(c(0, results[, 2:4])))
min.y.tick <- min(breaks_pretty()(c(0, results[, 2:4])))
my.ticks.y <- breaks_pretty()(c(0, results[, 2:4]))
label.stat1 <- stat1.name
label.stat2 <- stat2.name
label.cv <- paste(" cv(", 100 * alpha, "%)", sep = "")
if (scale.font != 1) {
scale.fontaxis <- 0.9 * scale.font
} else {
scale.fontaxis <- scale.font
}
scale.key <- 15 + 5 * scale.font
g.stat <- ggplot(df.results, aes(x = .data$lag)) +
theme_classic() +
theme(plot.title = element_text(size = 13 * scale.font, face = "bold"), legend.text = element_text(size = 12 * scale.font), axis.text = element_text(size = 10 * scale.fontaxis), axis.title = element_text(size = 10 * scale.font)) +
theme(legend.position = "top", legend.key.size = unit(scale.key, "pt"), legend.text = element_text(margin = margin(r = 10, unit = "pt"))) +
theme(axis.title.x = element_text(margin = margin(t = 10))) +
theme(axis.text.x = element_text(margin = margin(t = 4))) +
guides(fill = guide_legend(order = 1)) +
theme(axis.title.y = element_blank()) +
labs(x = "Lag") +
scale_x_continuous(breaks = my.ticks) +
scale_y_continuous(breaks = my.ticks.y, limits = c(min.y.tick, max.y.tick), expand = c(0, 0)) +
geom_blank(aes(y = max.y.tick)) +
geom_blank(aes(y = min.y.tick)) +
geom_line(aes(y = .data$stat1, colour = "stat1", linetype = "stat1"), linewidth = 1, na.rm = TRUE) +
geom_line(aes(y = .data$stat2, colour = "stat2", linetype = "stat2"), linewidth = 1, na.rm = TRUE) +
geom_line(aes(y = .data$cv, colour = "cv", linetype = "cv"), linewidth = 1, na.rm = TRUE)
if (stat.id == 1) {
g.stat <- g.stat + scale_colour_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "gray50", "stat2" = "#B32026", "cv" = "black"), labels = c(label.stat1, label.stat2, label.cv)) +
scale_linetype_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "solid", "stat2" = "solid", "cv" = "dashed"), labels = c(label.stat1, label.stat2, label.cv))
}
if (stat.id == 2 | stat.id == 3) {
g.stat <- g.stat + scale_colour_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "#B32026", "stat2" = "gray50", "cv" = "black"), labels = c(label.stat1, label.stat2, label.cv)) +
scale_linetype_manual("", breaks = c("stat1", "stat2", "cv"), values = c("stat1" = "solid", "stat2" = "solid", "cv" = "dashed"), labels = c(label.stat1, label.stat2, label.cv))
}
if (uni.biv == 1 && stat.id == 1) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cumulative tests for zero autocorrelation of ", my.names, sep = ""))
}
if (uni.biv == 1 && stat.id == 2) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Tests for iid property of ", my.names, sep = ""))
}
if (uni.biv == 1 && stat.id == 3) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cumulative tests for iid property of ", my.names, sep = ""))
}
if (uni.biv == 2 && stat.id == 1) {
g.stat <- g.stat + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(paste("Cumulative tests for zero cross-correlation of ", my.names[1], " and ", my.names[2], "(-k)", sep = ""))
}
show(g.stat)
}
myplotcorrmat <- function(cc, pv, my.names, scale.font, label.xi) {
nv <- NCOL(cc)
colnames(pv) <- my.names
rownames(pv) <- my.names
melted_pv <- melt(pv)
colnames(melted_pv) <- c("Var1", "Var2", "pv")
colnames(cc) <- my.names
rownames(cc) <- my.names
melted_cc <- melt(cc)
colnames(melted_cc) <- c("Var1", "Var2", "cc")
ccpv <- cbind(melted_cc, melted_pv[, 3])
colnames(ccpv) <- c("Var1", "Var2", "cc", "pv")
labelcc <- format(round(ccpv[, 3], 3), nsmall = 3)
labelcc[is.na(ccpv[, 4])] <- 1
labelpv <- paste("(", format(round(ccpv[, 4], 3), nsmall = 3), ")", sep = "")
labelpv[is.na(ccpv[, 4])] <- ""
mylabel <- paste(labelcc, "\n", labelpv, sep = "")
if (label.xi == TRUE) {
labelticks <- matrix(NA, nv, 1)
for (ii in 1:nv) {
labelticks[ii] <- as.expression(bquote(x[.(ii)]))
}
} else {
labelticks <- my.names
}
ccpv[is.na(ccpv)] <- 1.3
mycolors <- c("#B32026", "#C85C61", "#DC989B", "#F1D5D6", "#FFFFFF", "#BEBEBE")
colcode <- cut(ccpv[, 4], breaks = c(0, 0.001, 0.01, 0.05, 0.1, 1), include.lowest = TRUE)
brk <- levels(colcode)
colcode2 <- cut(ccpv[, 4], breaks = c(0, 0.001, 0.01, 0.05, 0.1, 1, 1.5), include.lowest = TRUE)
mylabelspv <- c("[0,0.001]", "(0.001,0.01]", "(0.01,0.05]", "(0.05,0.1]", "(0.1,1]")
scale.key <- 5 + 10 * scale.font
g <- ggplot(ccpv, aes(.data$Var1, rev(.data$Var2))) +
geom_tile(aes(fill = colcode2), colour = "black", show.legend = TRUE) +
scale_fill_manual(name = "p-value", breaks = brk, values = mycolors, labels = mylabelspv, drop = FALSE) +
scale_x_discrete(position = "top", labels = labelticks) +
scale_y_discrete(labels = rev(labelticks)) +
theme(plot.title = element_text(size = 13 * scale.font, face = "bold"), legend.title = element_text(size = 12 * scale.font), legend.text = element_text(size = 11 * scale.font), axis.text = element_text(size = 11 * scale.font)) +
theme(legend.key.size = unit(scale.key, "pt")) +
theme(plot.title = element_text(hjust = 0.5)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(panel.background = element_blank(), axis.ticks = element_blank()) +
theme(axis.ticks.length = unit(0, "pt")) +
theme(axis.text.x = element_text(colour = "black"), axis.text.y = element_text(colour = "black")) +
theme(plot.title = element_text(margin = margin(b = 20))) +
geom_text(aes(.data$Var1, rev(.data$Var2), label = mylabel), color = "black", size = 3.5 * scale.font) +
ggtitle("Pearson correlations and p-values for significance")
show(g)
}
Try the testcorr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.