R/inars_estimation.R

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

# This function is a adaptation of the betareg.control() function
# in betareg package
inars_control <- function(optim.method = "BFGS", maxit = 1000, hessian = TRUE,
                       trace = FALSE, start = NULL, ...)
{
  rval <- list(optim.method = optim.method, maxit = maxit, hessian = hessian,
               trace = trace, start = start)

  rval <- c(rval, list(...))

  if (!is.null(rval$fnscale))
    warning("fnscale must not be modified")

  rval$fnscale <- -1

  if (is.null(rval$reltol))
    rval$reltol <- .Machine$double.eps^(1/1.2)

  rval
}

## Initial values
inits <- function (innovation)
{
  switch(innovation,

         ## Skew discrete Laplace
         SDL = {
           start <- function(y){
             # Moment estimates
             alpha.h <- stats::acf(y, plot=F)$acf[2]
             c <- abs(alpha.h) * (1 - abs(alpha.h)) * mean(abs(y))
             mu.h <- (1 - alpha.h) * mean(y)

             if (1 - mu.h ^ 2 - 2 * c +
                 2 * (1 - alpha.h ^ 2) * stats::var(y) >= 0) {
               disp.h <- sqrt(1 - mu.h ^ 2 - 2 * c +
                                2 * (1 - alpha.h ^ 2) * stats::var(y)) - 1
             }else {
               warning("\nThe method of moments estimate of phi is not well defined for this sample\n")
               disp.h <- stats::var(y) * (1 - alpha.h ^ 2) + c + mu.h
             }

             c(alpha.h, mu.h, disp.h)
           }
         },

         ## Skellam
         SK = {
           start <- function(y){
             # Moment estimates
             alpha.h <- stats::acf(y, plot=F)$acf[2]
             c <- abs(alpha.h) * (1 - abs(alpha.h)) * mean(abs(y))
             mu.h <- (1 - alpha.h) * mean(y)
             disp.h <- (1 - alpha.h ^ 2) * stats::var(y) - c

             c(alpha.h, mu.h, disp.h)
           }
         },

         ## Discrete logistic
         DLOG = {
           start <- function(y){
             # Moment estimates
             alpha.h <- stats::acf(y, plot=F)$acf[2]
             c <- abs(alpha.h) * (1 - abs(alpha.h)) * mean(abs(y))
             mu.h <- (1 - alpha.h) * mean(y)
             q <- pi/sqrt(3 * ((1 - alpha.h ^ 2) * stats::var(y) - c - 0.0833))
             disp.h <- abs(exp(q)/(1 - exp(q)))

             c(alpha.h, mu.h, disp.h)
           }
         },

         stop(gettextf("%s innovation not recognised", sQuote(innovation)), domain = NA))

  environment(start) <- asNamespace("stats")
  structure(list(start = start, name = innovation))
}

# Parameter estimates of the INARS model
inars_est <- function(y, X = NULL, innovation = c("SDL", "SK", "DLOG"),
                      method = c("CML", "MM"),
                      optim.control = inars_control(...), ...)
{

  # Sample size
  n <- length(y)

  # Initial definitions
  innovation <- match.arg(innovation, c("SDL", "SK", "DLOG"))
  method <- match.arg(method, c("CML", "MM"))
  optim.method <- optim.control$optim.method
  maxit <- optim.control$maxit
  hessian <- optim.control$hessian
  trace <- optim.control$trace
  start  <- optim.control$start
  optim.control$method <- optim.control$hessian <- optim.control$start <- NULL

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

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

  if(is.null(start)){
    start <- inits(innovation)$start(y)
  }


  if (method == "MM") {
    if (p > 1L){
      stop("\nError: The moment estimator must be only used without regressors\n")
    }else{
      opt <- list()
      opt$par <- start
    }

  }else{

    # Conditional log--likelihood
    cll <- function(par){
      inars_cll(par, y, X, innovation = innovation)
    }

    # Conditional score
    if (innovation == "SDL"){
      cscore <- function(par){
        inars_cscore(par, y, X, innovation = innovation)
      }
    }else{
      cscore <- NULL
      optim.method <- "Nelder-Mead"
    }

    # Initial values
    if (p > 1) start <- c(start[1], solve(t(X)%*%X)%*%t(X)%*%y, start[3])


    opt <- suppressWarnings(stats::optim(par = start,
                                         fn = cll,
                                         gr = cscore,
                                         method = optim.method,
                                         control = optim.control,
                                         hessian = hessian))
    opt$start <- start

    if (opt$convergence > 0)
      warning(cat("optimization failed to converge\n"))
  }

  opt
}