R/sim.R

In PTSR: Positive Time Series Regression

#' @title
#' Function to simulate a PTSR model
#'
#' @details
#' The function \code{ptsr.sim} generates a random sample from a positive time
#' series regression model, with a given distribution.
#'
#' @param n a strictly positive integer. The sample size of yt (after burn-in).
#'  Default is 1.
#' @param burn a non-negative integer. length of "burn-in" period. Default is 0.
#' @param xreg optionally, a vector or matrix of external regressors.
#'  For simulation purposes, the length of xreg must be \code{n+burn}.
#'  Default is \code{NULL}.
#' @param xregar logical, if \code{FALSE}, the regressors are not included in the
#' AR component of the model. Default is \code{TRUE}.
#' @param varphi non-negative parameter. Default is 1.
#' @param alpha a numeric value corresponding to the intercept.  Default is 0.
#' @param beta optionally, a vector of coefficients corresponding to the
#'   regressors in \code{xreg}. Default is \code{NULL}.
#' @param phi optionally, for the simulation function this must be a vector
#'   of size \eqn{p}, corresponding to the autoregressive coefficients
#'   (including the ones that are zero), where \eqn{p} is the AR order.
#'   Default is \code{NULL}.
#' @param theta optionally, for the simulation function this must be a vector
#'   of size \eqn{q}, corresponding to the moving average coefficients
#'   (including the ones that are zero), where \eqn{q} is the MA order.
#'   Default is \code{NULL}.
#'  that \eqn{g_2(y_t) = 0}, for \eqn{t < 1}.
#' @param seed optionally, an integer which gives the value of the fixed
##'  seed to be used by the random number generator. If missing, a random integer
##'  is chosen uniformly from 1,000 to 10,000.
#' @param rdist function, the random number generator to be used
#' @param link1 character indicating which link must be used for \eqn{\mu_t}.  See \code{\link{ptsr.link}} for available links. Default is \sQuote{log}.
#' @param link2 character indicating which link must be used for \eqn{y_t} in the AR recursion.  See \code{\link{ptsr.link}} for available links. Default is \sQuote{identity}.
#' @param g1 optionally, a link function to be used  for \eqn{\mu_t}. Default is \code{NULL}, so that it is calculated internally, using \code{link1}.
#' @param g1.inv optionally, a the inverse link function to be used  for \eqn{\eta_t}. It must the the ivnerse of \code{g1}. Default is \code{NULL}, so that it is calculated internally, using \code{link1}.
#' @param g2 optionally, a link function to be used  for \eqn{y_t}. Default is \code{NULL}, so that it is calculated internally, using \code{link2}.
#'
#' @return
#' Returns a list with the following components
#' \itemize{
#'
#' \item \code{yt}: the simulated time series
#'
#' \item \code{mut}: the conditional mean
#'
#' \item \code{etat}: the linear predictor \eqn{g(\mu_t)}
#'
#' \item \code{error}: the error term.
#'
#' }
#'
#' @examples
#'
#' #-------------------------------------------------------------------
#' # Generating a sample of a Gamma-ARMA(1,1) model with no regressors
#' #-------------------------------------------------------------------
#'
#' simu = ptsr.sim(n = 300, burn = 50,
#'                 varphi = 20, alpha = 0,
#'                 phi = 0.35, theta = 0.2,
#'                 seed = 1234, rdist = r.gamma,
#'                 link1 = "log", link2 = "log")
#'
#' names(simu)
#' plot.ts(simu$yt)
#' lines(simu$mut, col= "red")
#'
#' @export
ptsr.sim <- function(n = 1, burn = 0, xreg = NULL, xregar = TRUE,
                     varphi = 1, alpha = 0, beta = NULL,
                     phi = NULL, theta = NULL,
                     seed = stats::runif(1, 1000, 10000), rdist = r.gamma,
                     link1 = "log", link2 = "identity",
                     g1 = NULL, g1.inv = NULL, g2 = NULL) {

  # function to be applied to mu and eta
  if (is.null(g1) | is.null(g1.inv)) {
    lk1 <- ptsr.link(link = link1)
    g1 <- lk1$linkfun
    g1.inv <- lk1$linkinv
  }

  # function to be applied to r2
  if (is.null(g2)) {
    lk2 <- ptsr.link(link = link2)
    g2 <- lk2$linkfun
  }

  set.seed(seed)

  p <- length(phi)
  q <- length(theta)

  n <- n + burn
  # vectors to save the time series
  yt <- numeric(n+p)
  gyt <- numeric(n+p)
  mut <- numeric(n)
  epst <- numeric(n+q)
  etat <- numeric(n)

  # checking beta and xreg
  if (is.null(beta) | is.null(xreg)) {
    XB <- matrix(0, nrow = n, ncol = 1)
  }
  else {
    if (is.null(xreg)) stop("missing xreg")
    xreg <- as.matrix(xreg)
    if (is.null(beta)) stop("missing beta")
    if (length(beta) != ncol(xreg)) stop("beta is not compatible with xreg")
    XB <- xreg %*% beta
  }
  if (nrow(XB) < n) stop("xreg must have n + burn rows")
  if(p > 0){
    # setting the initial p values as the unconditional mean
    if(xregar)	XB = c(rep(mean(XB[1:p]),p), XB)
  }

  # starting the loop
  for (t in 1:n) {

    etat[t] = alpha + XB[t]
    if(p > 0){
      iAR = (t+p) - (1:p)
      xr = 0
      if(xregar) xr = XB[iAR]
      etat[t] = etat[t] + sum(phi * (gyt[iAR] - xr))
    }
    if(q > 0) etat[t] = etat[t] + sum(theta * epst[(t+q) - (1:q)])

    mut[t] = g1.inv(etat[t])
    yt[t+p] = rdist(n = 1, mu = mut[t], varphi = varphi)
    gyt[t+p] = g2(yt[t+p])

    epst[t+q] = yt[t+p] - mut[t]
  }

	if(sum(yt[(p+burn + 1):(n+p)] <= 0) > 0) warning("Error generating the process. Zeros produced")

  invisible(list(yt = yt[(p+burn + 1):(n+p)],
                 mut = mut[(burn + 1):n],
                 etat = etat[(burn + 1):n],
                 epst = epst[(q+burn + 1):(n+q)]))
}