R/inars_cll.R

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

### Conditional log--likelihood function
inars_cll <- function(par, y, X = NULL,
                      innovation = c("SDL", "SK", "DLOG")){

  ## Design matrix
  if (is.null(X)) X <- matrix(1, nrow = length(y), ncol = 1)

  ## Sample size and dimension of the beta vector
  n <- length(y)
  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())

  ## Auxiliary function
  aux <- function(t){

    if (alpha * y[t-1] >= 0) {
      index <- (y[t] - abs(y[t-1])):(y[t])
      p <- log(sum(invt$d(index, mu[t], disp) *
                     stats::dbinom(y[t] - index, abs(y[t-1]), abs(alpha))))

    }else{
      index <- (y[t]):(abs(y[t-1]) + y[t])
      p <- log(sum(invt$d(index, mu[t], disp) *
                     stats::dbinom(index - y[t], abs(y[t-1]), abs(alpha))))

    }

    return(p)
  }

  # Conditional log--likelihood
  p <- apply(matrix(2:n, n - 1, 1), 1, aux)
  sum(p)
}

### Conditional score function
inars_cscore <- function(par, y, X = NULL,
                         innovation = c("SDL", "SK", "DLOG")){

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

  # Beta dimension vector and sample size
  p <- NCOL(X)
  n <- length(y)

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

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

  if (innovation != "SDL"){

    return(NULL)

  } else{

    invt <- match.fun(innovation)
    invt <- eval(invt())

    # Conditional probability function
    cpmf <- function(t){

      if (alpha * y[t-1] >= 0) {
        index <- (y[t] - abs(y[t-1])):(y[t])
        p <- sum(invt$d(index, mu[t], disp) *
                   stats::dbinom(y[t] - index, abs(y[t-1]), abs(alpha)))

      }else{
        index <- (y[t]):(abs(y[t-1]) + y[t])
        p <- sum(invt$d(index, mu[t], disp) *
                   stats::dbinom(index - y[t], abs(y[t-1]), abs(alpha)))

      }

      return(p)
    }

    ### For alpha (alpha != 0)
    dalpha <- function(k, t){
      choose(abs(y[t-1]), k) * alpha * (abs(alpha)^(k-2)) *
        ((1 - abs(alpha))^(abs(y[t-1]) - k - 1)) *
        (k - abs(y[t-1]) * abs(alpha))
    }

    ### For mu
    dmu <- function(k, t){
      id1 <- which(k >= 0); id2 <- which(k < 0)

      plus <- (((disp + mu[t])/(2 + disp + mu[t]))^(k[id1]-1)) *
        ((1/(2 + disp + mu[t]))^2) * 2 * k[id1]

      minus <- (((disp - mu[t])/(2 + disp - mu[t]))^(-k[id2]-1)) *
        ((1/(2 + disp - mu[t]))^2) * 2 * k[id2]

      dmu <- 1 / (1 + disp) * c(minus, plus)
      index2 <- c(which(k < 0), which(k >= 0))

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

    }


    ### For disp
    ddisp <- function(k, t){
      id1 <- which(k >= 0); id2 <- which(k < 0)

      minus <- -(((disp - mu[t])/(2 + disp - mu[t]))^(-k[id2]-1)) *
        ((1/(2 + disp - mu[t]))^2) * 2 * k[id2] - invt$d(k[id2], mu[t], disp)

      plus <- (((disp + mu[t])/(2 + disp + mu[t]))^(k[id1]-1)) *
        ((1/(2 + disp + mu[t]))^2) * 2 * k[id1] - invt$d(k[id1], mu[t], disp)

      ddisp <- c(minus, plus)/(disp + 1)
      index2 <- c(which(k < 0), which(k >= 0))

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

    ### Score function
    score <- function(t){
      if (alpha * y[t-1] >= 0) {
        index <- (y[t] - abs(y[t-1])):(y[t])
        ua <- sum(invt$d(index, mu[t], disp) *
                    dalpha(y[t] - index, t) / cpmf(t))
        ub <- sum(stats::dbinom(y[t] - index, abs(y[t-1]), abs(alpha)) *
                    dmu(index, t) / cpmf(t)) * X[t, ]
        up <- sum(stats::dbinom(y[t] - index, abs(y[t-1]), abs(alpha)) *
                    ddisp(index, t) / cpmf(t))
        c(ua, ub, up)

      }else{
        index <- (y[t]):(abs(y[t-1]) + y[t])
        ua <- sum(invt$d(index, mu[t], disp) *
                    dalpha(index - y[t], t) / cpmf(t))
        ub <- sum(stats::dbinom(index - y[t], abs(y[t-1]), abs(alpha)) *
                    dmu(index, t) / cpmf(t)) * X[t, ]
        up <- sum(stats::dbinom(index - y[t], abs(y[t-1]), abs(alpha)) *
                    ddisp(index, t) / cpmf(t))
        c(ua, ub, up)

      }

    }

    apply(apply(matrix(2:n, n - 1, 1), 1, score), 1, sum)

  }

}