R/LiMcLeod.R

In portes: Portmanteau Tests for Time Series Models

"LiMcLeod" <-
function(obj,lags=seq(5,30,5),fitdf=0,sqrd.res=FALSE){
  Maxlag <- max(lags)

  class.obj = class(obj)[1]
     TestType <- "0"
    if (class.obj == "ts" || class.obj == "numeric" || class.obj == "matrix" || class.obj == "mts")
        TestType <- "1"

    if (class.obj == "ar" || class.obj == "arima0" || class.obj ==
        "Arima" || class.obj == "ARIMA" || class.obj == "forecast_ARIMA" || class.obj == "varest" || class.obj == "list")
     TestType<-"2"

    if (TestType == "0")
        stop("obj must be class ar, arima0, Arima, (forecast_ARIMA ARIMA Arima), varest, ts, numeric, matrix, (mts ts), or list")

    if (TestType=="1")
       res <- as.ts(obj)
     else{
          GetResid <- GetResiduals(obj)
          res <- GetResid$res
          fitdf <- GetResid$fitdf
     }

     if (sqrd.res)
         res <- res^2
       n <- NROW(res)
       k <- NCOL(res)
	 df <- k^2*(lags-fitdf)
         NegativeDF <- which(df<0)
         df[NegativeDF] <- 0
	 Accmat <- stats::acf(res, lag.max = Maxlag, plot = FALSE, type = "correlation")$acf
	 inveseR0 <- solve(Accmat[1,,])
         prodvec <- numeric(Maxlag)
        for(l in 1:Maxlag){
           tvecR <- t(as.vector(Accmat[l+1,,]))
	   prodvec[l] <- crossprod(t(tvecR),crossprod(t(kronecker(inveseR0,inveseR0)),t(tvecR)))
	}
      Q <- n*cumsum(prodvec)
      STATISTIC <- Q[lags]+ k^2*lags*(lags+1)/(2*n)
      PVAL <- 1 - stats::pchisq(STATISTIC,df)
      PVAL[NegativeDF] <- NA
       summary <- matrix(c(lags,STATISTIC,df,PVAL),ncol=4)
       dimnames(summary) <- list(rep("", length(STATISTIC)),c("lags","statistic","df","p-value"))
    return(summary)
}