create_Z_varest.R

create_Z_varest <- function(object, ..., n.ahead = 1, ci = 0.95, dumvar = NULL) {
  
  # Extraction
  K <- object$K
  p <- object$p
  obs <- object$obs
  type <- object$type
  data.all <- object$datamat
  ynames <- colnames(object$y)
  n.ahead <- as.integer(n.ahead)
  Z <- object$datamat[, -c(1 : K)]
  B <- Bcoef(object)
  
  ##
  ## Deterministic and lagged y's
  ## Retrieval of A in matrix (whole)
  ## Deterministic variables in Zdet
  ##
  
  # Build trend component
  if(type == "const"){
    Zdet <- matrix(rep(1, n.ahead), nrow = n.ahead, ncol = 1)
    colnames(Zdet) <- "const"
  }else if(type == "trend"){
    trdstart <- nrow(Z) + 1 + p
    Zdet <- matrix(seq(trdstart, length = n.ahead), nrow = n.ahead, ncol = 1)
    colnames(Zdet) <- "trend"
  }else if(type == "both"){
    trdstart <- nrow(Z) + 1 + p
    Zdet <- matrix(c(rep(1, n.ahead), seq(trdstart, length = n.ahead)), nrow = n.ahead, ncol = 2)
    colnames(Zdet) <- c("const", "trend")
  }else if(type == "none"){
    Zdet <- NULL
  }
  
  ## Include seasonal if applicable
  if(!is.null(eval(object$call$season))){
    season <- eval(object$call$season)
    seas.names <- paste("sd", 1:(season-1), sep = "")
    cycle <- tail(data.all[, seas.names], season)
    seasonal <- as.matrix(cycle, nrow = season, ncol = season - 1)
    if(nrow(seasonal) >= n.ahead){
      seasonal <- as.matrix(cycle[1:n.ahead, ], nrow = n.ahead, ncol = season -1 )
    } else {
      while(nrow(seasonal) < n.ahead){
        seasonal <- rbind(seasonal, cycle)
      }
      seasonal <- seasonal[1:n.ahead, ]
    }
    rownames(seasonal) <- seq(nrow(data.all) + 1, length = n.ahead)
    if(!is.null(Zdet)){
      Zdet <- as.matrix(cbind(Zdet, seasonal))
    } else {
      Zdet <- as.matrix(seasonal)
    }      
  }
  
  ## Include exogenous variables if applicable
  if(!is.null(eval(object$call$exogen))){
    if(is.null(dumvar)){
      stop("\nNo matrix for dumvar supplied, but object varest contains exogenous variables.\n")
    }
    if(!all(colnames(dumvar) %in% colnames(data.all))){
      stop("\nColumn names of dumvar do not coincide with exogen.\n")
    }
    if(!identical(nrow(dumvar), n.ahead)){
      stop("\nRow number of dumvar is unequal to n.ahead.\n")
    }
    if(!is.null(Zdet)){
      Zdet <- as.matrix(cbind(Zdet, dumvar))
    } else {
      Zdet <- as.matrix(dumvar)
    }      
  }
  
  ## Retrieving predetermined y variables
  Zy <- as.matrix(object$datamat[, 1:(K * (p + 1))])
  
  ## forecast recursion
  forecast <- matrix(NA, ncol = K, nrow = n.ahead)
  lasty <- c(Zy[nrow(Zy), ])
  
  for(i in 1:n.ahead){
    lasty <- lasty[1 : (K * p)]
    Z <- c(lasty, Zdet[i, ])
    forecast[i, ] <- B %*% Z
    temp <- forecast[i, ]
    lasty <- c(temp, lasty)
  }
  
  Z <- as.data.frame(t(Z))
  names(Z) <- colnames(B)
  
  return(Z)
  
}