R/gFOBI.R

In tsBSS: Blind Source Separation and Supervised Dimension Reduction for Time Series

# Method gFOBI
gFOBI <- function(X, ...) UseMethod("gFOBI")

# main function for gFOBI
gFOBI.default <- function(X, k = 0:12, eps = 1e-06, maxiter = 100, method = c("frjd", "rjd"),
                           na.action = na.fail, weight = NULL,
                           ordered = FALSE, acfk = NULL, original = TRUE, alpha = 0.05, ...) {
  if (!is.numeric(X)) stop("non-numeric data")
  if (any(is.na(X) | is.infinite(X))) stop("missing/infinite values are not allowed")
  nk <- length(k)
  method <- match.arg(method)

  n <- nrow(X)
  p <- ncol(X)
  prep <- BSSprep(X)
  Y <- prep$Y 
  
  R <- array(0, dim = c(p, p, nk))
  for (i in 1:nk) {
    Yt <- Y[1:(n - k[i]), ]
    Yti <- Y[(1 + k[i]):n, ]
    r <- sqrt(rowSums(Yt^2))
    Yu <- r*Yti
    Ri <- crossprod(Yu)/nrow(Yt)
    R[, , i] <- Ri
  }
  JD <- switch(method, frjd = {
    JADE::frjd(R, eps = eps, maxiter = maxiter, na.action = na.action, weight = weight)$V
  }, rjd = {
    JADE::rjd(R, eps = eps, maxiter = maxiter, na.action = na.action)$V
  })
  W <- crossprod(JD, prep$COV.sqrt.i)
  W <- diag(sign(rowMeans(W))) %*% W
  S <- tcrossprod(prep$X.C, W)
  if (ordered == TRUE) { #Ordering by volatility
    if (is.null(acfk) == TRUE) { acfk <- k }
    ord <- ordf(S, acfk, p, W, alpha, ...)
    W <- ord$W
    if (original == TRUE) {
      S <- ord$S # Original independent components
    } else {
      S <- ord$RS # Residuals based on ARMA fit, if applicable; otherwise original IC's
      Sraw <- ord$S
      Sraw <- ts(Sraw, names = paste("Series", 1:p))
    }
  }
  S <- ts(S, names = paste("Series", 1:p))
  RES <- list(W = W, k = k, S = S, MU = prep$MEAN)
  if (ordered == TRUE) {
    if (original == FALSE) {
      RES$Sraw <- Sraw
    }
    RES$fits <- ord$fits
    RES$armaeff <- ord$armaeff
    RES$linTS <- ord$linTS
    RES$linP <- ord$linP
    RES$volTS <- ord$volTS
    RES$volP <- ord$volP
  }
  class(RES) <- c("bssvol", "bss")
  RES
}

gFOBI.ts <- function(X, ...) {
  x <- as.matrix(X)
  RES <- gFOBI.default(x, ...)
  S <- RES$S
  attr(S, "tsp") <- attr(X, "tsp")
  RES$S <- S
  if (!is.null(RES$Sraw)) {
    Sraw <- RES$Sraw
    attr(Sraw, "tsp") <- attr(X, "tsp")
    RES$Sraw <- Sraw
  }
  RES
}

gFOBI.xts <- function(X, ...) {
  x <- as.matrix(X)
  RES <- gFOBI.default(x, ...)
  S <- xts::as.xts(RES$S)
  attr(S, "index") <- attr(X, "index")
  xts::xtsAttributes(S) <- xts::xtsAttributes(X) #attributes additional to zoo
  RES$S <- S
  if (!is.null(RES$Sraw)) {
    Sraw <- xts::as.xts(RES$Sraw)
    attr(Sraw, "index") <- attr(X, "index")
    xts::xtsAttributes(Sraw) <- xts::xtsAttributes(X)
    RES$Sraw <- Sraw
  }
  RES
}

gFOBI.zoo <- function(X, ...) {
  x <- as.matrix(X)
  RES <- gFOBI.default(x, ...)
  S <- zoo::as.zoo(RES$S)
  attr(S, "index") <- attr(X, "index")
  RES$S <- S
  if (!is.null(RES$Sraw)) {
    Sraw <- zoo::as.zoo(RES$Sraw)
    attr(Sraw, "index") <- attr(X, "index")
    RES$Sraw <- Sraw
  }
  RES
}