R/leeg_simulation.R

In rdmatheus/leegarma:

#' LEEG-ARMA Generating of Bounded Time Series
#'
#' @description Generation of LEEG-ARMA process values with possible
#'  regression variables.
#'
#' @param n A strictly positive integer given the lenght of output series.
#' @param model A list with component \code{ar} and/or \code{ma} giving the
#'  AR and MA coefficients, respectively. Optionally, a component \code{beta}
#'  can be used to pass the coefficient vector associated with explanatory
#'  variables (without intercept).
#' @param X Optionally, a vector or matrix of explanatory variables,
#'  which must have number of rows equals to n.
#' @param alpha Intercept.
#' @param delta Dispersion parameter, which must be greater than zero.
#' @param freq The number of observations per unit of time (default: freq = 1).
#' @param link Character specification of the link function for the median
#'  model (mu). Currently, "logit", "probit", "cloglog", and "cauchit" are
#'  supported.
#'
#' @return Returns an object of class ts, which has an ARMA structure.
#' @export
#'
#' @author Diego R. Canterle
#' @author Rodrigo M. R. Medeiros <\email{rodrigo.matheus@live.com}>
leegarma_sim <- function(n, model = NULL, X = NULL,
                         alpha, delta,
                         freq = 1, link = "logit")
{
  if (!is.null(model) && !is.list(model))
    stop("'model' must be list, or must not be informed")
  if (n <= 0L)
    stop("'n' must be strictly positive")

  phi <- model$ar
  theta <- model$ma
  beta <- model$beta

  ar <- NULL
  ma <- NULL

  if(any(is.null(phi) == F))
  {
    ar <- 1:length(phi)
  }

  if(any(is.null(theta) == F))
  {
    ma <- 1:length(theta)
  }

  linktemp <- substitute(link)
  if (!is.character(linktemp))
  {
    linktemp <- deparse(linktemp)
    if (linktemp == "link")
      linktemp <- eval(link)
  }
  if (any(linktemp == c("logit", "probit", "cloglog","cauchit")))
  {
    stats <- make.link(linktemp)
  }else{
    stop(paste(linktemp, "link not available, available links are \"logit\", ",
               "\"probit\", \"cloglog\" and \"cauchit\""))
  }

  link <- structure(list(link = linktemp,
                         linkfun = stats$linkfun,
                         linkinv = stats$linkinv))

  linkfun <- link$linkfun
  linkinv <- link$linkinv


  ###### ARMA model
  if(any(is.null(phi) == F) && any(is.null(theta) == F) &&
         any(is.null(beta) == T))
  {
    p <- max(ar)
    q <- max(ma)
    m <- 2 * max(p,q)

    ynew <- rep(alpha, (n + m))
    mu <- linkinv(ynew)

    error <- rep(0, n + m) # E(error)=0
    eta <- y <- rep(NA, n + m)

    for(i in (m + 1):(n + m))
    {
      eta[i]  <- alpha + as.numeric(phi%*%ynew[i - ar]) +
                                  as.numeric(theta%*%error[i - ma])
      mu[i]   <- linkinv(eta[i])
      y[i]    <- rleeg(1, mu[i], delta)
      ynew[i] <- linkfun(y[i])
      error[i]<- ynew[i]-eta[i]
    }

    # ARMA model
    return( ts(y[(m + 1):(n + m)], frequency = freq) )
  }


  ###### AR model
  if(any(is.null(phi) == F) && any(is.null(theta) == T) &&
     any(is.null(beta) == T))
  {
    p <- max(ar)
    m <- 2 * p

    ynew <- rep(alpha, (n + m))
    mu <- linkinv(ynew)

    eta <- y <- rep(NA, n + m)

    for(i in (m + 1):(n + m))
    {
      eta[i]  <- alpha + (phi%*%ynew[i - ar])
      mu[i]   <- linkinv(eta[i])
      y[i]    <- rleeg(1, mu[i], delta)
      ynew[i] <- linkfun(y[i])
    }

    # AR model
    return( ts(y[(m + 1):(n + m)], frequency = freq) )
  }


  ###### MA model
  if(any(is.null(phi) == T) && any(is.null(theta) == F) &&
     any(is.null(beta) == T))
  {
    q <- max(ma)
    m <- 2 * q

    ynew <- rep(alpha, (n + m))
    mu <- linkinv(ynew)

    eta <- y <- error <- rep(0, n + m) # E(error)=0

    for(i in (m + 1):(n + m))
    {
      eta[i]   <- alpha + (theta%*%error[i - ma])
      mu[i]    <- linkinv(eta[i])
      y[i]     <- rleeg(1, mu[i], delta)
      ynew[i]  <- linkfun(y[i])
      error[i] <- ynew[i] - eta[i]
    }

    # MA model
    return( ts(y[(m+1):(n+m)],frequency=freq) )
  }

  ###### ARMAX model
  if(any(is.null(phi) == F) && any(is.null(theta) == F) &&
     any(is.null(beta) == F))
  {
    p <- max(ar)
    q <- max(ma)
    m <- 2 * max(p, q)

    ynew <- rep(alpha, (n + m))
    mu <- linkinv(ynew)

    error <- rep(0, n + m) # E(error)=0
    eta <- y <- rep(NA, n + m)

    X <- as.matrix(X)
    X <- rbind(matrix(0, ncol = dim(X)[2], nrow = m), X)

    for(i in (m+1):(n+m))
    {
      eta[i]   <- alpha + X[i,]%*%as.matrix(beta) +
        (phi%*%(ynew[i - ar] - X[i - ar,]%*%as.matrix(beta))) +
        (theta%*%error[i - ma])
      mu[i]    <- linkinv(eta[i])
      y[i]     <- rleeg(1, mu[i], delta)
      ynew[i]  <- linkfun(y[i])
      error[i] <- ynew[i] - eta[i]
    }

    # ARMAX model
    return(ts(y[(m + 1):(n + m)], frequency = freq))
  }

  ###### ARX model
  if(any(is.null(phi) == F) && any(is.null(theta) == T) &&
     any(is.null(beta) == F))
  {
    p <- max(ar)
    m <- 2 * p

    ynew <- rep(alpha, (n + m))
    mu <- linkinv(ynew)

    eta <- y <- rep(NA, n + m)

    X <- as.matrix(X)
    X <- rbind(matrix(0, ncol = dim(X)[2], nrow = m), X)

    for(i in (m+1):(n+m))
    {
      eta[i]  <- alpha + X[i,]%*%as.matrix(beta) +
        (phi%*%(ynew[i - ar] - X[i - ar, ]%*%as.matrix(beta)))
      mu[i]   <- linkinv(eta[i])
      y[i]    <- rleeg(1, mu[i], delta)
      ynew[i] <- linkfun(y[i])
    }
    # ARX model
    return(ts(y[(m + 1):(n + m)], frequency = freq))
  }

  ###### MAX model
  if(any(is.null(phi) == T) && any(is.null(theta) == F) &&
     any(is.null(beta) == F))
  {
    q <- max(ma)
    m <- 2 * q

    ynew <- rep(alpha, (n + m))
    mu <- linkinv(ynew)

    error <- rep(0, n + m) # E(error)=0
    eta<- y <- rep(NA, n + m)

    X <- as.matrix(X)
    X <- rbind(matrix(0, ncol = dim(X)[2], nrow = m), X)

    for(i in (m+1):(n+m))
    {
      eta[i]   <- alpha + X[i,]%*%as.matrix(beta) + (theta%*%error[i - ma])
      mu[i]    <- linkinv(eta[i])
      y[i]     <- rleeg(1, mu[i], delta)
      ynew[i]  <- linkfun(y[i])
      error[i] <- ynew[i] - eta[i]
    }

    # MAX model
    return(ts(y[(m+1):(n+m)],frequency=freq))
  }

}