R/XgenSimple.R

In yymmhaha/PackPaper1: Package related to my first paper

#' Generate multiple time series given simple VAR structure
#'
#' @description Generate vector autoregression time series with user-specified correlation structure
#'
#' @details This function has been used as the most basic data generator. It is also used in the \code{XgenCorr, XgenCorr2}
#' and intervention functions.
#' Mathematical details see the book.
#'
#' @param p a number. Number of covariates x.
#' @param A a list. Each element is a transition (lag coefficient) matrix. Caution: some specifications can make the series unstationary
#' @param x0 a list. Each element is a vector of length p, initial value of x. The length depends on the time lags (number of A).
#' @param t a number. Length of series, including the initial lags.
#' @param sigmax a number. Standard deviation of noise of x.
#'
#' @return a list of components
#'
#' \item{Big A}{The expanded transition matrix for all lags}
#' \item{X}{The generated multiple time series }
#' \item{lag}{The lag for VAR process}
#' \item{Characteristic roots}{Absolute value for 1/eigen(BigA), if any is smaller than 1, process is unstationary and series can not be generated}
#'
#'
#' @export
#'
#' @examples
#' A <- list(A1 = toeplitz(c(0.4, rep(0, 9))),
#'          A2 = toeplitz(c(0.2, rep(0, 9))),
#'          A3 = toeplitz(c(0.1, rep(0, 9))))
#' x0 <- list(x1 = rep(0, 10), x2 = rep(0, 10), x3 = rep(0, 10))
#' sigmax <- 1
#' nseries <- 10
#' nsize <- 50
#'
#' XgenSimple(p = nseries, A = A, x0 = x0, t = nsize, sigmax = sigmax)


XgenSimple <- function(p, A, x0, t, sigmax){

  lagx <- length(A)
  d <- lagx * p                   # dimension of the big matrix

  # A part
  BigA <- matrix(rep(0, d^2), d, d)
  A.top <- do.call(cbind, A)
  if(lagx == 1){                 # when there's only 1 lag
    BigA <- A.top
  }else{
    D <- diag(1, p*(lagx-1))
    BigA[1:p, ] <- A.top
    BigA[(p+1):(p*lagx), 1:(p*(lagx-1))] <- D
  }

  # check stationarity
  charoot <- abs(1/eigen(BigA)$values)
  if(length(charoot[charoot<=1])!=0){
    message(paste("Process not stationary. Change your coefficient matrices"))
    return(list('BigA' = BigA, 'Characteristic roots' = charoot))
  }else{

    # X part
    BigX <- matrix(rep(0, p*t*lagx), (p*lagx), t)
    X0 <- do.call(c, x0)
    BigX[, 1] <- X0

    for (i in 1:(t-1)){
      noise <- c(rnorm(p, 0, sigmax), rep(0, (p*(lagx-1))))
      BigX[, (i+1)] <- BigA %*% BigX[, i] + noise
    }
    X.final <- t(BigX[1:p, ])  # take the top p rows and transpose!

    return(list('BigA' = BigA, 'X' = X.final, 'lag' = lagx,
                'Characteristic roots' = charoot))
  }
}