R/inars_innovations.R

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

#' @name innovations
#'
#' @title Family of Innovations Process for Fitting a INARS(1) model
#'
#' @description Provide the current available distributions that can be
#'     used as an innovation process in a fit of the INARS(1) model.
#'
#' @return The above functions return a list of the following components:
#' \describe{
#'   \item{d}{ Probability mass function \code{function(z, mu, disp)}.}
#'   \item{r}{ A random generator \code{function(n, mu, disp)}.}
#' }
#'   The parameters associated with these functions are:
#' \itemize{
#'   \item{\code{z}:}{ An integer-valued quantile.}
#'   \item{\code{mu}:}{ The real-valued mean of the process innovation.}
#'   \item{\code{disp}:}{ The positive dispersion parameter of the process innovation.}
#'   \item{\code{n}:}{ Numbers of random values to return.}
#' }
#'
#' These functions are used to define the innovation process in the \code{inars()} fit.
#'
#'
#'
#'
#' @author Rodrigo M. R. Medeiros <\email{rodrigo.matheus@live.com}>
#'
#' @references
#' Skellam, J. G. (1946). The frequency distribution of the difference
#'     between two Poisson variates belonging to different populations.
#'     Journal of the Royal Statistical Society: Series A, 109, 296.
#'
#' Andersson, J., & Karlis, D. (2014). A parametric time series model with
#'     covariates for integers in Z. Statistical Modelling, 14, 135--156.
#'
#' Chakraborty, S., & Chakravarty, D. (2016). A new discrete probability
#'     distribution with integer support on (--infty, infty). Communications in
#'     Statistics-Theory and Methods, 45, 492--505.
#'
#' Medeiros, R.M.R., Bourguignon M. (2021) A useful regression model for
#'     integers on Z with application to quality of life in
#'     prisons. Submitted.
#'
#'
#' @examples
#' sdl <- SDL()
#' skellam <- SK()
#' dlog <- DLOG()
#'
#' zvals <- -10:10
#'
#' ## SDL
#' barplot(sdl$d(zvals, -2, 4), names.arg = zvals)
#' barplot(sdl$d(zvals, 0, 4), names.arg = zvals)
#' barplot(sdl$d(zvals, 2, 4), names.arg = zvals)
#'
#' ## Skellam
#' barplot(skellam$d(zvals, -2, 4), names.arg = zvals)
#' barplot(skellam$d(zvals, 0, 4), names.arg = zvals)
#' barplot(skellam$d(zvals, 2, 4), names.arg = zvals)
#'
#' ## Discret logistic
#' barplot(dlog$d(zvals, -2, 1), names.arg = zvals)
#' barplot(dlog$d(zvals, 0, 1), names.arg = zvals)
#' barplot(dlog$d(zvals, 2, 1), names.arg = zvals)
#'
NULL

### Skew discrete Laplace innovation process
#' @rdname innovations
#' @export
SDL <- function(){
  out <- list()

  ## Probability mass function
  out$d <- function(z, mu, disp){

    pmf <- rep(0, length(z))
    if (length(mu) == 1)  mu <- rep(mu, length(z))
    if (length(disp) == 1)  disp <- rep(disp, length(z))

    # NaN index
    NaNindex <- which(abs(mu) >= disp | disp <= 0)
    pmf[NaNindex] <- NaN

    # Positive pmf index
    index <- which(z %% 1 == 0 & !is.nan(pmf))

    prob.minus <- ((disp[index & z < 0] - mu[index & z < 0]) /
                     (2 + disp[index & z < 0] - mu[index & z < 0])) ^
      abs(z[which(z %% 1 == 0 & !is.nan(pmf) & z < 0)])
    prob.plus  <- ((mu[index & z >= 0] + disp[index & z >= 0]) /
                     (2 + mu[index & z >= 0] + disp[index & z >= 0])) ^
      z[which(z %% 1 == 0 & !is.nan(pmf) & z >= 0)]

    pmf[index] <- 1 / (1 + disp[index]) * c(prob.minus, prob.plus)
    index2 <- c(which(z < 0), which(z >= 0))

    return(pmf[sort(index2, index.return=T)$ix])
  }

  ## Random generation
  out$r <- function(n, mu, disp){
    if (disp <= 0)
      warning("\nThe dispersion parameter must be positive.\n")
    if (any(disp < abs(mu)))
      warning("\nConstraints are not satisfied.\n")

    return(stats::rgeom(n, 2 / (2 + disp + mu)) -
           stats::rgeom(n, 2 / (2 + disp - mu)))
  }

  out
}

### Skellam innovation process
#' @rdname innovations
#' @export
SK <- function(){
  out <- list()

  ## Probability mass function
  out$d <- function(z, mu, disp)
  {
    pmf <- rep(0, length(z))
    if (length(mu) == 1)  mu <- rep(mu, length(z))
    if (length(disp) == 1)  disp <- rep(disp, length(z))

    # NaN index
    NaNindex <- which(abs(mu) >= disp | disp <= 0)
    pmf[NaNindex] <- NaN

    # Positive pmf index
    index <- which(z %% 1 == 0 & !is.nan(pmf))

    theta1 <- (disp[index] + mu[index])/2
    theta2 <- (disp[index] - mu[index])/2
    exp(-(theta1 + theta2)) * (theta1/ theta2) ^ (z[index]/2) *
      besselI(2 * sqrt(theta1 * theta2), nu = abs(z[index]))
  }

  # Random generation
  out$r <- function(n, mu, disp){
    if (disp <= 0)
      warning("\nThe standard deviation must be positive.\n")
    if (any(disp < abs(mu)))
      warning("\nConstraints are not satisfied.\n")

    return(stats::rpois(n, (disp + mu)/2) -
           stats::rpois(n, (disp - mu)/2))
  }

  out
}

### Discret logistic innovation
#' @rdname innovations
#' @export
DLOG <- function(){
  out <- list()

  ## Probability mass function
  out$d <- function(z, mu, disp)
  {
    pmf <- rep(0, length(z))
    if (length(mu) == 1)  mu <- rep(mu, length(z))
    if (length(disp) == 1)  disp <- rep(disp, length(z))

    # NaN index
    NaNindex <- which(disp <= 0)
    pmf[NaNindex] <- NaN

    # Positive pmf index
    index <- which(z %% 1 == 0 & !is.nan(pmf))

    mu <- mu + 0.5
    p <- disp / (1+disp)

    (1 - p) * (p^(z - mu)) /
    ((1 + p^(z - mu)) * (1 + p^(z - mu + 1)))
  }

  ## Random generation
  out$r <- function(n, mu, disp){
    if (disp <= 0)
      warning("\nThe standard deviation must be positive.\n")

    mu <- mu + 0.5
    p <- disp / (1 + disp)

    if ((n == 1L) & (length(mu) > 1)) mu <- mu[1]

    u <- stats::runif(n)
    ceiling(log((1-u)/u)/log(p) - 1 + mu)
  }

  out
}