R/decor.R

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

##' The function decorrelates a Data Set so that the new data will have no correlation or very small correlations
##'
##' The input x must be a vector. decor() decorrelate a Data Set based on Cholesky decomposition. It is assumed that the data observations are covariance stationary and the serial correlation exists only when two observations are within bmax > 0 in their observation indices. More specifically, it is assumed that γ(q) = Cov(Xi,Xi+q) only depends on q when i changes, and γ(q) = 0 when q > bmax, where Xi and Xi+q are two process observations obtained at times i and i + q. In practice, the autocorrelation between Xi and Xi+q usually decays when q increases. In such cases, γ(q) is small when q is large, and thus a proper value of bmax can be chosen such that γ(q) ≈ 0 when q > bmax. The default value of bmax is 10, which indicates that the assumption is γ(q) = Cov(Xi,Xi+q) = 0 when q is greater than 10
##' @title Decorrelate data observations
##' @param x univariate data
##' @param bmax The smallest number that any data points have no or small correlations with previous data points which means γ(q) ≈ 0 when q > bmax
##' @return The data after decorrelation
##' @author Xiulin Xie
##' @export
##' @import matrixcalc
##' @examples
##' decor(rnorm(100,0,1),bmax = 9)

decor <- function(x, bmax = 10) {
    cov_vector <- Vcov(x, bmax)
    cov_matrix <- Vcovmatrix(x, bmax)
    avg <- mean(x)
    nx <- rep(0, length(x))
    for (i in 1:length(x)) {
        if (i == 1) {
            nx[i] <- (x[i] - avg)/cov_vector[i]^0.5
            b = 1
        } else {
            if (i == 2) {
                d <- cov_vector[1] - cov_vector[2]/cov_matrix[1, 1] * cov_vector[2]
                e <- x[1] - avg
                nx[i] <- (x[i] - avg - cov_vector[2]/(cov_matrix[1, 1]) * e)/sqrt(d)
                b <- min(b + 1, bmax)
            } else {
                d <- cov_vector[1] - cov_vector[(b + 1):2] %*% matrix.inverse(cov_matrix[1:b, 1:b]) %*% cov_vector[(b + 
                  1):2]
                e <- x[(i - b):(i - 1)] - avg
                nx[i] <- (x[i] - avg - cov_vector[(b + 1):2] %*% matrix.inverse(cov_matrix[1:b, 1:b]) %*% 
                  e)/sqrt(d)
                b <- min(b + 1, bmax)
            }
        }
        
    }
    return(nx)
}