R/MVcovmatrix.R

In XiulinXie/SPCmonitor2: Process Control Chart for Monitoring Short-Range Serially Correlated Data

##' The function MVcovmatrix computes the Complete Covariance Matrix of a multivariate data set x and the covariance or correlation matrix between x(i) and x(i+b). Covariance γ(q) = Cov(Xi,Xi+q) of a Data Set where X is a multivariate data set
##'
##' The input x must be a matrix. MVcovmatrix() computes complete covariance matrix between X(i) and X(i+q),where X(i) and X(i+q) are multivariate process observations obtained at times i and i + q where q =0,1,2.....bmax. Covariance γ(q) = Cov(Xi,Xi+q), where X is a p-dimensinal vector and γ(q) will be a p by p matrix. The default value of bmax is 10, which means the output will be [[γ(0),γ(1),.....γ(10)],.....[γ(i-1),γ(i),.....γ(11-i),...[γ(10),γ(9),.....γ(0)]] where γ(q) = Cov(X(i),X(i+q)) a p by p matrix, so the output will be a 11*p by 11*p matrix.
##' @title The Complete Covariance Matrix of a Multivariate Data Set
##' @param x data
##' @param bmax γ(q) = Cov(Xi,Xi+q) , the maximun value of q
##' @return  Multivariate Covariance Matrix M.
##' @author Xiulin Xie
##' @export
##' @examples
##' MVcovmatrix(matrix(rnorm(900,0,1),nrow = 3))
MVcovmatrix <- function(x, bmax = 10) {
    cov_vector <- MVcov(x)$cov_vector
    p <- nrow(x)
    cov_matrix = matrix(rep(0, p * p * (bmax + 1) * (bmax + 1)), nrow = nrow(x) * (bmax + 1))
    for (i in 1:(bmax + 1)) {
        for (j in 1:(bmax + 1 - i + 1)) {
            l = ((i - 1) * nrow(x) + 1)
            o1 = ((j - 1) * nrow(x) + 1)
            o2 = ((j + i - 2) * nrow(x) + 1)
            cov_matrix[o1:(o1 + p - 1), o2:(o2 + p - 1)] = cov_vector[, l:(l + p - 1)]
            cov_matrix[o2:(o2 + p - 1), o1:(o1 + p - 1)] = t(cov_vector[, l:(l + p - 1)])
        }
    }
    return(cov_matrix)
}