R/inars_simulation.R

In rdmatheus/inars: INARS(1) Model for Integer-Valued Time Series

#' @name sim-inars
#'
#' @title INARS(1) Process Random Generation
#'
#' @param n number of random values to return.
#' @param par parameter vector; it must be specified in the following order:
#'  (alpha, mu, disp) to simulate without regressors; or (alpha, beta, disp)
#'  to simulate with a regression structure in the mean of the innovation
#'  process. In the latter, beta represents a \code{p} dimensional vector of
#'  the coefficients.
#' @param X design matrix for the mean. If NULL (default), the parameter
#'  vector \code{par} must be specified as (alpha, mu, disp).
#' @param innovation the assumed distribution for the innovation process.
#'  Currently, are available the skew discrete Laplace ("SDL"), the Skellam
#'  ("SK"), and the discret logistic ("DLOG") distributions.
#'
#' @return A integer-values time series of size n, which consists of a
#' realization of the INARS(1) process with an innovation process specified
#' via \code{innovation} argument.
#' @export
#'
#' @references
#' Kim, H. Y., & Park, Y. (2008). A non-stationary integer-valued
#'     autoregressive model. Statistical papers, 49, 485.
#'
#' Andersson, J., & Karlis, D. (2014). A parametric time series model with
#'     covariates for integers in Z. Statistical Modelling, 14, 135--156.
#'
#' @author Rodrigo M. R. Medeiros <\email{rodrigo.matheus@live.com}>
#'
#' @examples
#' \dontrun{
#' # Sample size
#' n <- 100
#'
#' ###############################
#' # Generate without regressors #
#' ##############################
#'
#' # Parameters
#' alpha <- -0.5
#' mu <- -3
#' disp <- 4
#'
#' # SDL innovations
#' y <- inars_sim(n, c(alpha, mu, disp))
#' barplot(table(y), xlab = "y", ylab = "Frequency")
#' plot(y, xlab = "Time", ylab = "y")
#'
#' # SK innovation
#' y <- inars_sim(n, c(alpha, mu, disp), innovation = "SK")
#' barplot(table(y), xlab = "y", ylab = "Frequency")
#' plot(y, xlab = "Time", ylab = "y")
#'
#' ################################################
#' # Generate with a regression structure in the  #
#' #      mean of the innovation process          #
#' ###############################################
#'
#' # Parameters
#' alpha <- -0.5
#' beta <- c(1.2, 2)
#' disp <- 4
#'
#' # Design matrix
#' X <- cbind(rep(1, n), runif(n, -1, 1))
#'
#' # SDL innovations
#' y <- inars_sim(n, c(alpha, beta, disp), X)
#' barplot(table(y), xlab = "y", ylab = "Frequency")
#' plot(y, xlab = "Time", ylab = "y")
#'
#' # SK innovations
#' y <- inars_sim(n, c(alpha, beta, disp), X, innovation = "SK")
#' barplot(table(y), xlab = "y", ylab = "Frequency")
#' plot(y, xlab = "Time", ylab = "y")
#' }
#'
inars_sim <- function(n, par, X = NULL,
                      innovation = c("SDL", "SK", "DLOG")){

  if (is.null(X)) X <- matrix(1, nrow = n, ncol = 1)

  # Dimension of the beta vector
  p <- NCOL(X)

  # Parameters
  alpha <- par[1]
  mu <- X%*%as.matrix(par[2:(p + 1)])
  disp <- par[p + 2]

  # Innovation process distribution
  invt <- match.fun(match.arg(innovation, c("SDL", "SK", "DLOG")))
  invt <- eval(invt())

  y <- vector()
  nn <- 100 + n
  if (n < 100){
    mu <- c(rep(mu, ceiling(nn/n))[1:100], mu)
  } else {
    mu <- c(mu[1:100], mu)
  }

  y[1] <- invt$r(1, mu[1], disp)
  z <- invt$r(nn, mu, disp)

  # Signed binomial thinning operation
  sign_thinning <- function(alpha, y){
    sign(alpha) * sign(y) * stats::rbinom(1, abs(y), abs(alpha))
  }

  # Random generation
  for (t in 2:nn){
    y[t] <- sign_thinning(alpha, y[t-1]) + z[t]
  }

  return(stats::ts(y[101:nn]))
}