R/weighted_ordinal_pattern_distribution.R

In statcomp: Statistical Complexity and Information Measures for Time Series Analysis

## Sebastian Sippel
# 05.01.2014


#' @title A function to compute weighted ordinal pattern statistics
#' @export
#' @useDynLib statcomp, .registration = TRUE
#' @description Computation of weighted ordinal patterns of a time series. Weights can be generated by a user-specified function (e.g. variance-weighted, see Fadlallah et al 2013).
#' @usage weighted_ordinal_pattern_distribution(x, ndemb)
#' @param x A numeric vector (e.g. a time series), from which the weighted ordinal pattern distribution is to be calculated  
#' @param ndemb Embedding dimension of the ordinal patterns (i.e. sliding window size). Should be chosen such as length(x) >> ndemb
#' @details 
#' This function returns the distribution of weighted ordinal patterns using the Keller coding scheme, detailed in Physica A 356 (2005) 114-120. NA values are allowed.
#' The function uses old and slow R routines and is only maintained for comparability. 
#' For faster routines, see \link{weighted_ordinal_pattern_distribution}.
#' @return A character vector of length factorial(ndemb) is returned.
#' @references Fadlallah, B., Chen, B., Keil, A. and Principe, J., 2013. Weighted-permutation entropy: A complexity measure for time series incorporating amplitude information. Physical Review E, 87(2), p.022911.
#' @author Sebastian Sippel
#' @seealso \code{\link{weighted_ordinal_pattern_distribution}}
#' @examples
#' x = arima.sim(model=list(ar = 0.3), n = 10^4)
#' weighted_ordinal_pattern_distribution(x = x, ndemb = 6)
weighted_ordinal_pattern_distribution = function(x, ndemb) {
  
  
  epsilon=1.e-10
  npdim=factorial(ndemb)
  
  #Berechnungs der Ordnungsstatistik nach der Kodierung von Karsten Keller:
  #Physica A 356 (2005) 114-120
  
  nfac=factorial(ndemb)
  
  wifrec=.C("weighted_ordinal_pattern_loop",
            as.double(x),
            as.integer(length(x)),
            as.integer(ndemb),
            numeric(nfac),
            as.integer(nfac), NAOK=T)[[4]]
    
  # ifrec is the ordinal pattern distribution in the Keller coding scheme!
  return(wifrec)
}


#' @keywords internal
weighted_ordinal_pattern_distribution_2 = function(x, ndemb, weight.fun = var.fun) {
  
  ### Deal with gaps in the sliding window time series:
  # get indices to run through for calculation of complexity measures
  gapfree = sapply(1:(length(x)-ndemb + 1), FUN = function(y) if(!any(is.na(x[y:(y+ndemb-1)]))) return(y))
  
  epsilon=1.e-10
  npdim=factorial(ndemb)
  
  #Berechnungs der Ordnungsstatistik nach der Kodierung von Karsten Keller:
  #Physica A 356 (2005) 114-120
  
  ifrec = numeric(length=npdim)  #ersetzt die for-schleife zum erstellen des Vektors, #for ip=1:npdim; ifrec( ip ) = 0; end;
  w_ifrec = numeric(length=npdim) # Vektor f?r weigthed Permutations: see Fadlallah et al (2013), Eq. 4
  N_weighted = numeric(length=1) # Variable that stores weights, see Fadlallah et al (2013), Eq. 4
  
  ## introduce for loop:
  for (nv in 1:(length(gapfree))) {
    
    xvec <- x[gapfree[nv]:(gapfree[nv] + ndemb - 1)]
    
    ## only if no gaps are in the "word" of the time series:
    ipa = matrix(data=0,nrow=ndemb, ncol=ndemb)  #Inversionsmatrix
    
    for (il in 2:ndemb) {
      for (it in il:ndemb) { 
        ipa[it, il] = ipa[it-1, il-1]
        if( (xvec[it] <= xvec[it - ( il - 1 ) ] ) || ( abs( xvec[it - ( il - 1)] - xvec[it]) < epsilon))
          ipa[ it, il ] = ipa[ it, il ] + 1;
      }
    }
    
    nd = ipa[ndemb,2] 
    for (il in 3:ndemb) {
      nd =il * nd + ipa[ndemb, il]
    }
    
    ifrec[nd + 1] = ifrec[nd+ 1] + 1;  # ordinal pattern distribution (non-weighted!)     
    
    # calculation of weights for weighted permutation entropy:
    weights_xvec <- weight.fun(xvec)
    w_ifrec[nd + 1] = w_ifrec[nd + 1] + 1 * weights_xvec #Adds weights to current permutation pattern
    # N_weighted = N_weighted + 1*weights_xvec #Calculates denominator of Eq. 4 in Fadlallah et al (2013)
    # N_weighted is simply the sum of w_ifrec!
  }
  
  # ifrec is the ordinal pattern distribution in the Keller coding scheme!
  return(w_ifrec)
}