#' 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))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.