R/autocovariances.R

Defines functions rgarch1p1 nvcovOfAcfBD nvcovOfAcf acfOfSquaredArmaModel ciplot plot_acr cilines plot.acf.test whiteNoiseTest acfGarchTest wnacf_asyse wnacf_asycov_garch wnacf_asycov_iid acfMaTest acfIidTest comboAcf comboAcvf backwardPartialVariances partialVariances partialAutocovariances partialAutocorrelations autocovariances basecfclass2S4 printFittedPart

Documented in acfGarchTest acfIidTest acfMaTest acfOfSquaredArmaModel autocovariances backwardPartialVariances nvcovOfAcf nvcovOfAcfBD partialAutocorrelations partialAutocovariances partialVariances rgarch1p1 whiteNoiseTest

## Do not edit this file manually.
## It has been automatically generated from *.org sources.

setClass("VirtualAutocovariances", contains = "VIRTUAL")
setClass("VirtualAutocorrelations", contains = "VIRTUAL")

setClass("Autocovariances", contains = c("FlexibleLagged", "VirtualAutocovariances"))
setClass("PartialAutocovariances", contains = c("FlexibleLagged", "VirtualAutocovariances"))
setClass("Autocorrelations", contains = c("FlexibleLagged", "VirtualAutocovariances"))
setClass("PartialAutocorrelations", contains = c("FlexibleLagged", "VirtualAutocovariances"))

## Commenting out the slots, their intention is valid in the univariate case
setClass("PartialVariances", contains = c("FlexibleLagged", "VirtualAutocovariances")
         ## , slots = c(data = "numeric")
        )

## Combo variants
setClass("ComboAutocovariances", contains = c("FlexibleLagged", "VirtualAutocovariances")
         ## , slots = c(data = "matrix")
        )

setClass("ComboAutocorrelations", contains = c("FlexibleLagged", "VirtualAutocorrelations")
         ## , slots = c(data = "matrix")
        )

setClass("Fitted",
         slots = c(n = "numeric", varnames = "character", objectname = "character" )
         )

setClass("SampleAutocovariances",         contains = c("Autocovariances",        "Fitted"))
setClass("SamplePartialAutocovariances",  contains = c("PartialAutocovariances", "Fitted"))
setClass("SampleAutocorrelations",        contains = c("Autocorrelations",       "Fitted"))
setClass("SamplePartialAutocorrelations", contains = c("PartialAutocorrelations","Fitted"))
setClass("SamplePartialVariances",        contains = c("PartialVariances",       "Fitted"))
## TODO: FittetComboXX variants?

.printFittedPart <- function(x){
    for(slot in slotNames("Fitted")){
        cat('Slot ', slot, ':', "\n", sep = "")
        print(slot(x, slot))
    }
}

setMethod("show", signature(object = "Autocovariances"),
          function (object)
          {
              .reportClassName(object, "Autocovariances")
              callNextMethod()
          }
          )

setMethod("show", signature(object = "Autocorrelations"),
          function (object)
          {
              .reportClassName(object, "Autocorrelations")
              callNextMethod()
          }
          )

setMethod("show", signature(object = "SampleAutocorrelations"),
          function (object)
          {
              .reportClassName(object, "SampleAutocorrelations")
              callNextMethod()
              .printFittedPart(object)
          }
          )

setMethod("show", signature(object = "SampleAutocovariances"),
          function (object)
          {
              .reportClassName(object, "SampleAutocovariances")
              callNextMethod()
              .printFittedPart(object)
          }
          )

setMethod("show", signature(object = "SamplePartialVariances"),
          function (object)
          {
              .reportClassName(object, "SamplePartialVariances")
              callNextMethod()
              .printFittedPart(object)
          }
          )

.basecfclass2S4 <- function(x, r0){
    stopifnot(class(x) == "acf")

    lagged <- Lagged(x) # Lagged knows how to interpret "acf" objects
    if(!missing(r0))
        lagged[0] <- if(class(r0) == "acf") r0$acf[1, , ] else r0

    data <- lagged[]
    n <- x$n.used
    series <- x$series
    varnames <- if(is.null(x$snames)) character(0)  else x$snames #TODO: NA in the former?
    ## lag <- x$lag

    class <- switch(x$type, covariance  = "SampleAutocovariances",
                            correlation = "SampleAutocorrelations",
			    partial     = "SamplePartialAutocorrelations" )

    new(class, data = data, n = n, varnames = varnames, objectname = series)
}

autocovariances <- function(x, maxlag, ...){
    if(is.list(x) && all(names(x) %in% c("ar", "ma", "sigma2"))){ # theoretical ARMA autocov
                 # ARMAacf(x$ar, x$ma, lag.max = maxlag)
                 #    TODO: if using ARMAacf, need to calculate R(0) and multiply by it!
                 #          using the function from package ltsa instead - but need checking;
                 #          does it work for non-invertible models?
        ## TODO: currently only univariate

        phi    <- if(is.null(x$ar     )) numeric(0) else x$ar
        theta  <- if(is.null(x$ma     )) numeric(0) else x$ma
        sigma2 <- if(is.null(x$sigma2)) NA_real_   else x$sigma2

                                    # !!! Note the '-' below!
        wrk <- tacvfARMA(phi = phi, theta = - theta, sigma2 = sigma2, maxLag = maxlag)

        new("Autocovariances", data = wrk)
    }else{ ## numeric, "ts", "mts" etc. ...
                           # TODO: check what happens if argument 'plot' is  in `...' (error?)
        wrk <- if(missing(maxlag))
                   acf(x, plot = FALSE, type = "covariance", ...)
               else
                   acf(x, lag.max = maxlag, plot = FALSE, type = "covariance", ...)
        obj <- .basecfclass2S4(wrk)

                  # TODO: FUNCTIONS IN PCTS PACKAGE WHICH RELY ON RESULT OF TYPE slMatrix()
                  #       IN THE MULTIVARIATE CASE;
                  #       THEY SHOULD CONTERT TO IT THEMSELVES.
                  #
                  # res <- acfbase2sl(res)
                  # slMatrix(res)
        obj
    }
}

setGeneric("autocovariances")

setGeneric("autocorrelations",
           function(x, maxlag, lag_0, ...)  standardGeneric("autocorrelations")
           )

partialAutocorrelations <- function(x, maxlag, lag_0 = TRUE, ...){
            # TODO: needs more work!
            #       see handwritten notes "parcor"; James Proberts may wish to do it.
            #
            # in this case, compute autocovariances, solve YW and extract the parcor.
            #   (for the multivariate case needs modification!)
            #
            # calculate autocor not autocov, since the latter may need more, e.g. sigma2

    acrobj <- if(missing(maxlag))
                  autocorrelations(x, ...)
              else
                  autocorrelations(x, maxlag = maxlag, ...)
    pacr <- modelCoef(acrobj, "PartialAutocorrelations")

    if(isTRUE(lag_0)){
        if(is(acrobj, "SampleAutocorrelations")){
             new("SamplePartialAutocorrelations", data = pacr, as(acrobj, "Fitted"))
        }else
             new("PartialAutocorrelations", data = pacr)
    }else if(identical(lag_0, FALSE)){
            # if(dim([email protected]) == 3)
            #     new("FlexibleLagged", data = [email protected][ -1, , ])
            # else
            #     new("FlexibleLagged", data = [email protected][-1])
        ## TODO: this is very lazy; to avoid checking the dimensions
        ##       do it properly (maybe via an auxiliary function)
        wrk <- new("FlexibleLagged", data = pacr)
        wrk[1:maxLag(wrk)]
    }else {# lag_0 = "var" or anything else - put the variance in lag_0
            ## TODO: this may fail if 'x' doesn't contain info about autocovariances.
            ##    Maybe just put NA instead?
            wrk <- autocovariances(x, maxlag = 0, ...)
            res <- as(obj, "FlexibleLagged")
            res[[0]] <- wrk[[0]]
            res
    }

}

setGeneric("partialAutocorrelations", signature = c("x", "maxlag", "lag_0") )

partialAutocovariances <- function(x, maxlag, ...){
            # TODO: see handwritten notes "parcor"!
            # TODO: for now don't use lag_0 since it is not supported consistently yet.
            # wrk <- autocorrelations(x, lag_0 = "var", ...) # TODO: put R(0) at lag(0)
        # TODO: only univariate for now
    acvfobj <- autocovariances(x, maxlag, ...)
    pacvf <- .comboAcvf(acvfobj, "pacvf")

    if(is(acvfobj, "SampleAutocovariances")){
        new("SamplePartialAutocovariances", data = pacvf, as(acvfobj, "Fitted"))
    }else
        new("PartialAutocovariances", data = pacvf)
}

setGeneric("partialAutocovariances", signature = c("x", "maxlag") )

partialVariances <- function(x, ...){ # improvement? was: "predictionVariances"
                        # TODO: see handwritten notes "parcor"!
                        #       see also comments in the default for partialAutocovariances()
         # TODO: only univariate for now
    acvfobj <- autocovariances(x, ...)
    wrk <- .comboAcvf(acvfobj, "psigma2")

    new("PartialVariances", data = wrk)
    if(is(acvfobj, "SampleAutocovariances")){
        new("SamplePartialVariances", data = wrk, as(acvfobj, "Fitted"))
    }else
        new("PartialVariances", data = wrk)
}

setGeneric("partialVariances") # was: setGeneric("predictionVariances")

backwardPartialVariances <- function(x,...) stop(x)
setGeneric("backwardPartialVariances")      # setGeneric("backwardPredictionVariances")

partialCoefficients <-
backwardPartialCoefficients <- function(x, p) stop(x)
setGeneric("partialCoefficients")	       # setGeneric("predictionCoefficients")	       
setGeneric("backwardPartialCoefficients")   # setGeneric("backwardPredictionCoefficients")

## TODO: set default value for maxlag?
setMethod("autocovariances",
    signature(x = "VirtualArmaModel"),
    function (x, maxlag, ...)
    {
        co <- modelCoef(x, convention = "BJ", ...)
        sigma2 <- sigmaSq(x) # method for this how to call it:
                             #      sigmaSq(), innovationVariance(), innovationVar()?
        if(missing(maxlag))
            maxlag <- max(length(co$ar), length(co$ma))
                ## TODO: compare to
                ##       FitARMA::TacvfARMA(phi = co$ar, theta = co$theta, lag.max = maxlag)
        res <- ltsa::tacvfARMA(phi = co$ar, theta = co$ma, maxLag = maxlag,
                               sigma2 = sigma2)
        as(res, "Autocovariances")
    }
)

setMethod("autocorrelations",
    signature(x = "ANY", maxlag = "ANY", lag_0 = "ANY"),
    function (x, maxlag, lag_0, ...)
    {
        obj <- autocorrelations(x, maxlag, ...)

        if(isTRUE(lag_0))
            obj
        else if(identical(lag_0, FALSE)){
            ## drop lag 0 and the class
	        # if(!is.null(dim([email protected])) && dim([email protected]) == 3)
                #     [email protected][ -1, , ]
	        # else
                #     [email protected][-1]
            obj[seq(length = maxLag(obj))] # was:  obj[][-1]
        }else {# lag_0 = "var" or anything else - put the variance in lag_0
                ## TODO: this may fail if 'x' doesn't contain info about autocovariances.
                ##    Maybe just put NA instead?
                wrk <- autocovariances(x, maxlag = 0, ...)
                res <- as(obj, "FlexibleLagged")
                res[[0]] <- wrk[[0]]
                res
        }

    }
)

setMethod("autocorrelations",
          signature(x = "ANY", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              if(isS4(x))
                  stop("there is no applicable method for this object's class")

              if(is.list(x) && all(names(x) %in% c("ar", "ma", "sigma2"))){
                  ## theoretical ARMA autocov
                  wrk <- ARMAacf(x$ar, x$ma, lag.max = maxlag, ...) # `...' here may be only
                                                                    # lag.max TODO: take care
                                                                    # of lag_0
                  new("Autocorrelations", data = wrk)
              }else{ # numeric, "ts", "mts" etc. ...
                  ## sample autocorrelations
                  wrk <- if(missing(maxlag))
                      acf(x, plot = FALSE, type = "correlation", ...)
                  else
                      acf(x, lag.max = maxlag, plot = FALSE, type = "correlation", ...)
                  .basecfclass2S4(wrk)
              }
              ## TODO: RESULT CHANGED FROM THAT IN pcts, ACCOMODATE FOR THIS IN pcts!
              ##
              ## res <- acfbase2sl(res) # convert to season/lag format
              ## res <- slMatrix(res)
              ## if(isTRUE(lag_0))
              ##     res
              ## else if(identical(lag_0, FALSE))
              ##     as.matrix(res)[-1]

              ## else {# lag_0 = "var" or anything else - put the variance in lag_0
              ##     stop('lag_0 = "var" not implemented yet for this method.')
              ## }
          }
          )



## seemingly redundant method but it does some non-trivial work ('maxlag' and possibly changes
## the class of x.
setMethod("autocorrelations",
          signature(x = "Autocorrelations", maxlag = "missing", lag_0 = "missing"),
          function (x, ...)
          {
              as(x, "Autocorrelations")   # x may be from a subclass
          }
          )

setMethod("autocorrelations",
          signature(x = "Autocorrelations", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, lag_0, ...)
          {
              new("Autocorrelations", data = x[0:maxlag]) # may introduce NA's
          }
          )


setMethod("autocorrelations",
          signature(x = "Autocovariances", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              acr <- x[0:maxlag] / x[[0]] # TODO: Only univariate case!
              new("Autocorrelations", data = acr)
          }
          )

setMethod("autocorrelations",
          signature(x = "PartialAutocovariances", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              acr <- modelCoef(x, "Autocorrelations")
              res <- new("Autocorrelations", data = acr)
              if(!missing(maxlag))
                  res@data <- res[0:maxlag]     # may introduce NA's
              res
          }
          )

setMethod("autocorrelations",
          signature(x = "PartialAutocorrelations", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              acr <- modelCoef(x, "Autocorrelations")
              res <- new("Autocorrelations", data = acr)
              if(!missing(maxlag))
                  res@data <- res[0:maxlag]     # may introduce NA's
              res
          }
          )


setMethod("autocorrelations",
          signature(x = "SamplePartialAutocovariances", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              acr <- modelCoef(x, "Autocorrelations")
              res <- new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
              if(!missing(maxlag))
                  res@data <- res[0:maxlag]     # may introduce NA's
              res
          }
          )

setMethod("autocorrelations",
          signature(x = "SamplePartialAutocorrelations", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              acr <- modelCoef(x, "Autocorrelations")
              res <- new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
              if(!missing(maxlag))
                  res@data <- res[0:maxlag]     # may introduce NA's
              res
          }
          )

## TODO: set default value for maxlag?
setMethod("autocorrelations",
    signature(x = "VirtualArmaModel", maxlag = "ANY", lag_0 = "missing"),
    function (x, maxlag, ...)
    {
        ## lazy - calls the default method with "BD" convention
        co <- modelCoef(x, convention = "BD", ...)
        co$sigma2 <- 1 # for autocorrelations this doesn't matter, but may be NA in x
        if(missing(maxlag))
            maxlag <- max(length(co$ar), length(co$ma))

        autocorrelations(co, maxlag = maxlag, ...)
    }
)

setMethod("autocorrelations",
    signature(x = "VirtualSarimaModel", maxlag = "ANY", lag_0 = "missing"),
    function (x, maxlag, ...)
    {
        arma <- as(x, "ArmaModel") # should give error if x is non-sttionary
        autocorrelations(arma, maxlag = maxlag, ...)
    }
)

setMethod("partialAutocorrelations",
          signature(x = "ts", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              wrk <- if(missing(maxlag))
                  pacf(x, plot = FALSE)
              else
                  pacf(x, lag.max = maxlag, plot = FALSE)
              .basecfclass2S4(wrk)
          }
          )

setMethod("partialAutocorrelations",
          signature(x = "mts", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...){             # TODO: check what definition they use, etc.
              wrk <- if(missing(maxlag))
                  pacf(x, plot = FALSE)
              else
                  pacf(x, lag.max = maxlag, plot = FALSE)
              r0 <- acf(x, lag.max = 0, plot = FALSE)
              res <- .basecfclass2S4(wrk, r0)
              if(!missing(maxlag))
                  res@data <- res[0:maxlag]     # may introduce NA's
          }
          )


setMethod("partialAutocorrelations",
          signature(x = "PartialAutocovariances", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              pacr <- modelCoef(x, "PartialAutocorrelations")
              res <- new("PartialAutocorrelations", data = pacr)
              if(!missing(maxlag))
                  res@data <- res[0:maxlag]     # may introduce NA's
              res
          }
          )

## not needed since the default does the job
##
## setMethod("partialAutocorrelations",
##     signature(x = "VirtualArmaModel"),
##     function (x, maxlag, lag_0 = TRUE, ...)
##     {
##         co <- armaCoef(x, convention = "BJ", ...)
##         sigma2 <- sigmaSq(x)
##
##         wrk <- FitARMA::TacvfARMA(phi = numeric(0), theta = numeric(0), lag.max = maxlag)
##         ## TODO: mult by sigma2
##     }
## )

.comboAcvf <- function(acvf, ind){ # acvf is a Lagged object (actually Autocovariance or similar)
    acvf <- acvf[] # drop the "Laggedness" (i.e. get the raw data)

    ## scalar case only; for now
    r0 <- acvf[1]

    acf <- acvf/r0
    acfonly <- as(acf, "vector")[-1]

    wrk <- ltsa::DLAcfToAR(acfonly) # further args: (useC = TRUE, PDSequenceTestQ = FALSE)
          ## wrk is a matrix with columns containing ar, pacf, sigma2;
              ##   (but sigma2 (column "sigsqk") contains standardised innov vars,
              ##    i.e. ones corresponding to R[0]=1);
          ## modify wrk below to our requirements.

          ## innovation variances in wrk are for R[0] = 1, so multiple them by R[0]
    wrk[ , 3] <- r0 * wrk[ , 3]
    wrk[ , 2] <- wrk[ , 2] * wrk[ , 3] # pacf => pacvf

    res <- rbind(c(1, r0, r0), wrk) # prepend a row  (1 R[0] R[0])'
         # 2016-11-24 was: res <- cbind([email protected], res)           # prepend also the acvf's
    res <- cbind(acvf[], res)           # prepend also the acvf's
                           # 'p' in 'pacvf' and 'psigma2' is for 'partial' as in pacf, etc.
    colnames(res) <- c("acvf", "ar", "pacvf", "psigma2")
    rownames(res) <- 0 : maxLag(acvf)

    res <- t(res)   # 2016-12-07 now return the lags as columns!

    if(missing(ind))
        res
    else
        res[ind, ]
}

## not used?
.comboAcvfNames <- c("Autocovariances"        = "acvf",
                     "ar"                     = "ar", # not used
                     "PartialAutocovariances" = "pacvf",
                     "PartialVariances"       = "sigma2" )


.comboAcf <- function(acvf, ind){ # acvf or acf
    acvf <- acvf[] # drop the "Laggedness" (i.e. get the raw data)

    ## scalar case only; for now
    r0 <- acvf[1]

    acf <- acvf/r0
    acfonly <- acf[-1]

    wrk <- ltsa::DLAcfToAR(acfonly) # further args: (useC = TRUE, PDSequenceTestQ = FALSE)
          ## wrk is a matrix with columns containing ar, pacf, sigma2;
              ##   (but sigma2 (column "sigsqk") contains standardised innov vars,
              ##    i.e. ones corresponding to R[0]=1);
          ## modify wrk below to our requirements.

    res <- cbind(acfonly, wrk)       # prepend also the acf's
    res <- rbind(c(1, 1, 1, 1), res) # prepend a row of 1's
    colnames(res) <- c("acf", "ar", "pacf", "stdsigma2")
    rownames(res) <- 0 : maxLag(acvf)

    res <- t(res)   # 2016-12-07 now return the lags as columns!

    if(missing(ind))
        res
    else
        res[ind, ]
}


## setAs("Autocovariances", "ComboAutocovariances",
##       function(from){
##
##
##           new("ComboAutocovariances", coef = )
##       }
## )

setAs("vector", "Autocovariances", function(from) new("Autocovariances", data = from))
setAs("vector", "Autocorrelations", function(from) new("Autocorrelations", data = from))
setAs("vector", "PartialAutocovariances",
      function(from) new("PartialAutocovariances", data = from))
setAs("vector", "PartialAutocorrelations",
      function(from) new("PartialAutocorrelations", data = from))

setAs("ANY", "PartialAutocovariances",
      function(from){
          partialAutocovariances(from)
      })

setAs("ANY", "PartialVariances",
      function(from){
          partialVariances(from)
      })

setAs("ANY", "Autocorrelations",
      function(from){
          autocorrelations(from)
      })

setAs("ANY", "PartialAutocorrelations",
      function(from){
          partialAutocorrelations(from)
      })

setAs("ANY", "ComboAutocovariances",
      function(from){
          acvf <- as(from, "Autocovariances")
          wrk <- .comboAcvf(acvf)
          new("ComboAutocovariances", data = wrk)
      })

setAs("ANY", "ComboAutocorrelations",
      function(from){
          acf <- as(from, "Autocorrelations")
          wrk <- .comboAcf(acf)
          new("ComboAutococorrelations", data = wrk)
      })

setAs("Autocovariances", "ComboAutocovariances",
      function(from){
          wrk <- .comboAcvf(from)
          new("ComboAutocovariances", data = wrk)
      })

setAs("Autocovariances", "ComboAutocorrelations",
      function(from){
          wrk <- .comboAcf(from)
          new("ComboAutocorrelations", data = wrk)
      })

setAs("Autocorrelations", "ComboAutocovariances",
      function(from){
          stop("not possible, missing R(0)")
      })

setAs("Autocorrelations", "ComboAutocorrelations",
      function(from){
          wrk <- .comboAcf(from)
          new("ComboAutocorrelations", data = wrk)
      })

setAs("PartialVariances", "Autocovariances",
      function(from){
          stop("PartialVariances cannot be converted to the requested class.")
      })

setAs("PartialVariances", "ComboAutocovariances",
      function(from){
          stop("PartialVariances cannot be converted to the requested class.")
      })

setAs("PartialVariances", "Autocorrelations",
      function(from){
          stop("PartialVariances cannot be converted to the requested class.")
      })

setAs("PartialVariances", "ComboAutocorrelations",
      function(from){
          stop("PartialVariances cannot be converted to the requested class.")
      })

setMethod("modelCoef", c("VirtualAutocovariances", "missing", "missing"),
          function(object){
              object@data
          }
          )

setMethod("modelCoef", c("VirtualAutocovariances", "character", "missing"),
          function(object, convention){
              if(identical(class(object), convention))
                  ## this case is equivalent to the one-argument call()
                  modelCoef(object)
              else{
                  convention <- new(convention)
                  modelCoef(object, convention  = convention )
              }
          }
          )

setMethod("modelCoef", c("VirtualAutocovariances", "VirtualAutocovariances", "missing"),
          function(object, convention){
              acvf <- modelCoef(object, "Autocovariances")
              modelCoef(acvf, convention)
          }
          )

setMethod("modelCoef", c("Autocovariances", "ComboAutocovariances", "missing"),
          function(object, convention){
              .comboAcvf(object)
          }
          )

setMethod("modelCoef", c("Autocovariances", "ComboAutocorrelations", "missing"),
          function(object, convention){
              .comboAcf(object)
          }
          )

# not defined since R[0] is unknown
# setMethod("modelCoef", c("Autocorrelations", "ComboAutocovariances", "missing") )

setMethod("modelCoef", c("Autocorrelations", "ComboAutocorrelations", "missing"),
          function(object, convention){
              .comboAcf(object)
          }
          )

setMethod("modelCoef", c("ComboAutocovariances", "Autocovariances", "missing"),
          function(object, convention){
              object@data["acvf", ]
          }
          )

## TODO: to AR model (or filter) using $ar

setMethod("modelCoef", c("ComboAutocovariances", "PartialAutocovariances", "missing"),
          function(object, convention){
              object@data["pacvf", ]
          }
          )

setMethod("modelCoef", c("ComboAutocovariances", "PartialVariances", "missing"),
          function(object, convention){
              object@data["sigma2", ]
          }
          )


## acvf to acf etc.
setMethod("modelCoef", c("ComboAutocovariances", "VirtualAutocovariances", "missing"),
          function(object, convention){
              acvf <- modelCoef(object, "Autocovariances")
              modelCoef(acvf, convention)
          }
          )

setMethod("modelCoef", c("ComboAutocorrelations", "Autocorrelations", "missing"),
          function(object, convention){
              object@data["acf", ]
          }
          )

## TODO: to AR model (or filter) using $ar

setMethod("modelCoef", c("ComboAutocorrelations", "PartialAutocorrelations", "missing"),
          function(object, convention){
              object@data["pacf", ]
          }
          )

## need class StandardizedPartialVariances
## setMethod("modelCoef", c("ComboAutocorrelations", "StandardizedPartialVariances", "missing"),
##           function(object, convention){
##           [email protected]["stdsigma2", ]
##       }
##       )

setMethod("modelCoef", c("Autocovariances", "Autocorrelations", "missing"),
          function(object, convention){
              co <- object[] # or modelCoef(object) - no need for such generality?
              co / co[1]        # TODO: modify to work in multivariate case
          }
          )

setMethod("modelCoef", c("Autocovariances", "PartialAutocorrelations", "missing"),
          function(object, convention){
                    # alternatively, obtain them from combo (which uses ltsa::DLAcfToAR;
                    # using here acf2AR; compare to combo during testing.
              co <- object[] # or modelCoef(object) - no need for such generality?
              ar.all <- acf2AR(co) # all partial prediction coef's (YW)
              c(1, diag(ar.all) )
          }
          )

## same as for "Autocovariances"
setMethod("modelCoef", c("Autocorrelations", "PartialAutocorrelations", "missing"),
          function(object, convention){
              co <- object@data # or modelCoef(object) - no need for such generality?
              ar.all <- acf2AR(co) # all partial prediction coef's (YW)
              c(1, diag(ar.all) )
          }
          )

setMethod("modelCoef", c("PartialAutocorrelations", "Autocorrelations", "missing"),
          function(object, convention){
                  # 2016-11-24 was:
                  #     pacf <- [email protected]
                  #     pacfonly <- pacf[-1]
              pacfonly <- object[1:maxLag(object)]
              ar <- FitAR::PacfToAR(pacfonly)
                   # alternatively:
                   # acvf <- ltsa::tacvfARMA(phi = ar, maxLag = length(ar))
                   # acf <- acvf/acvf[0]
              ARMAacf(ar, lag.max = length(ar))
          }
          )

##not available on purpose
##setMethod("modelCoef", c("Autocorrelations", "PartialAutocovariances", "missing"))

setMethod("modelCoef", c("PartialAutocovariances", "PartialAutocorrelations", "missing"),
          function(object, convention){
              pacr <- object[] / object[0] # TODO: currently scalar only
          }
          )

## n     - length of time series, scalar
## npar  - no. estimated param, scalar
## nlags - number of acf for the LB test, vector of positive integers
##
acfIidTest <- function(acf, n, npar = 0, nlags = npar + 1,
                       method = c("LiMcLeod", "LjungBox", "BoxPierce"),
                       interval = 0.95, expandCI = TRUE, ...){
    method <- match.arg(method)
    maxlag <- max(nlags)
    usedLags <- 1:maxlag
    switch(method,
           "LjungBox" = {
               rsq <- acf[usedLags]^2 / (n - usedLags)
               Q <- n * (n+2) * cumsum(rsq)
           },
           "BoxPierce" = {
               rsq <- acf[usedLags]^2
               Q <- n * cumsum(rsq)
           },
           "LiMcLeod" = {
               rsq <- acf[usedLags]^2
               Q <- n * cumsum(rsq) + usedLags * (usedLags + 1) / (2*n)
           },
           stop("Unknown method")
           )

    Q <- Q[nlags] # keep only the requested ChiSq values

    df <- nlags - npar
    df[df <= 0] <- 1    # TODO: this should really be NA or depend on the estimation method

    pval <- pchisq(Q, df = df, lower.tail = FALSE)

    wrk <- cbind(ChiSq = Q, DF = df, pvalue = pval)
    attr(wrk, "method") <- method
    res <- list(test = wrk)

    if(!is.null(interval)){
        lq <- qnorm((1 - interval) / 2)
        int <- lq / sqrt(n)
        if(expandCI)
            int <- rep(int, maxlag)

        int <- cbind(int, - int)
        res$ci <- int
	attr(res$ci, "level") <- interval
    }

    res
}

setGeneric("acfIidTest")

setMethod("acfIidTest", "SampleAutocorrelations",
          function(acf, n, npar, nlags, method, interval = 0.95, x){
              if(missing(n))
                  n <- acf@n
              else
                  if(n != acf@n){
                      ## only warning to allow 'what-if use.
                      warning("argument 'n' is not equal to slot 'n' of the autocorrelations")
                  }
              callNextMethod(acf, n = n, nlags = nlags, method = method, interval = interval)
          })

setMethod("acfIidTest", "missing",
          function(acf, n, npar, nlags, method, interval = 0.95, x){
              if(missing(n))
                  n <- length(x)
              acf <- autocorrelations(x, maxlag = max(nlags))
              acfIidTest(acf, n = n, nlags = nlags, method = method, interval = interval)
          })

acfMaTest <- function(acf, ma,  n, nlags,
                       interval = 0.95){
    maxlag <- max(nlags)
    usedLags <- 1:maxlag

    W <- nvcovOfAcfBD(acf, ma = ma, maxlag = maxlag)

    Q <- numeric(length(nlags))
    for(i in seq(along = nlags)){
        h <- nlags[i]
        la <- (ma + 1):h
        rho <- acf[la]
        stat <- rho %*% solve(W[la,la]) %*% rho
        Q[i] <- stat / n
    }

    df <- nlags - ma

    pval <- pchisq(Q, df = df, lower.tail = FALSE)

    wrk <- cbind(ChiSq = Q, DF = df, pvalue = pval)
    attr(wrk, "method") <- paste0("Ma(", ma, ") Test")
    res <- list(test = wrk)

    if(!is.null(interval)){
        lq <- qnorm((1 - interval) / 2, sd = sqrt(diag(W)))
        int <- lq / sqrt(n)

        int <- cbind(int, - int)
        res$ci <- int
        attr(res$ci, "level") <- interval
    }

    res
}

## 'x' is the time series
.wnacf_asycov_iid <- function(x, maxlag, forcor = TRUE, n = length(x)){
    if(forcor)
        diag(1/n, nrow = maxlag)
    else{
        r0 <- acf(x, type = "covariance", lag.max = 0, plot = FALSE)$acf[1, , ]
        diag(r0^2/length(x), nrow = maxlag) # TODO: verify!
    }
}

.wnacf_asycov_garch <- function(x, maxlag, forcor = TRUE){ # TODO: default for maxlag
    x <- as.vector(x) # otherwise x*wrk  below doesn't work
    n <- length(x)
    wrk <- sapply(1:maxlag, function(i) c(rep(0, i), x[1:(n-i)]))
    wrk <- x * wrk # multiple each column by x using recycling
    res <- acf(wrk, type = "covariance", lag.max = 0, plot = FALSE, demean = FALSE)
    res <- res$acf[1, , ] # extract the zero lag matrix
    if(forcor){
        r0 <- acf(x, type = "covariance", lag.max = 0, plot = FALSE,
                  demean = FALSE)$acf[1, , ]
        res <- res / r0^2
    }
    res / n  # since sqrt(n)* rhat has asy.cov. 'res'
}

.wnacf_asyse <- function(x, maxlag, h0 = "iid", forcor = TRUE, n = length(x)){
    switch(h0,
           "iid" = {
               v <- .wnacf_asycov_iid(x = x, maxlag = maxlag, forcor = forcor, n = n)
           },
           "garch" = {
               v <- .wnacf_asycov_garch(x = x, maxlag = maxlag, forcor = forcor)
           },
           ##default
           stop("Unrecognised null hypothesis")
           )
    ses <- sqrt(diag(v))
    lags <- 1:maxlag

    cbind(Lag = lags, StdError = ses) #  todo: include 'Estimate = cf[parm]'
}

acfGarchTest <- function(acr, x, nlags, interval = 0.95){
    v <- .wnacf_asycov_garch(x, maxlag = max(nlags))
         # the formulas in the book multiply by n, since Gamma is 'n times asy.cov.'
         #    n <- length(x)
         #    Q <- n * sapply(nlags, function(h) acr[1:h] %*% solve(v[1:h,1:h]) %*% acr[1:h] )

    Q <- sapply(nlags, function(h) acr[1:h] %*% solve(v[1:h,1:h]) %*% acr[1:h] )

    pval <- pchisq(Q, df = nlags, lower.tail = FALSE)
    res <- list(test = cbind(h = nlags, Q = Q, pval = pval))

    if(!is.null(interval)){
        lq <- qnorm((1 - interval) / 2)
        se <- sqrt(diag(v))
        int <- lq * se
        int <- cbind(int, -int)
        res$ci <- int
    }
    res
}

whiteNoiseTest <- function(object, h0, ...){
    switch(h0,
           "iid" = {
               acfIidTest(object, ...)
           },
           "garch" = {
               acfGarchTest(object, ...) # no method yet for this case
           },
           ##default
           stop("Unrecognised null hypothesis")
           )
}

setMethod("coef", "SampleAutocorrelations",
          function (object, ...){
              dataWithLagNames(object)
          }
          )

setMethod("vcov", "SampleAutocorrelations",
          function (object, assuming = "iid", maxlag = maxLag(object), ...) {
              res <- switch(assuming,
                            iid = .wnacf_asycov_iid(n = object@n, maxlag = maxlag),
                            garch = .wnacf_asycov_garch(maxlag = maxlag, ...),
                            stop("Unknown assumption")
                            )
              if(is.null(names(res)))
                  ## TODO: currently scalar case only;
                  colnames(res) <- rownames(res) <- paste0("Lag_", 1:nrow(res))
              res
          }
          )

# setMethod("diagOfVcov", "SampleAutocorrelations",

setMethod("confint", "SampleAutocorrelations",
          function (object, parm, level = 0.95, se = FALSE, maxlag, ...){
              ## based on confint.default(), main difference: we pass "..." to vcov();
              ## arg. 'se'
              cf <- coef(object)[-1] # drop lag 0
              pnames <- names(cf)
              if (missing(parm)){
                  if(missing(maxlag)){
                      lags <- 1:length(pnames)
                      parm <- pnames
                  }else{
                      lags <- 1:maxlag
                      parm <- pnames[1:maxlag]
                  }
              }else if (is.numeric(parm)){
                  lags <- parm
                  parm <- pnames[parm]
              }
              a <- (1 - level)/2
              a <- c(a, 1 - a)
              pct <- .stats.format.perc(a, 3)
              fac <- qnorm(a)
              ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, pct))
                   # TODO: pass maxlag to vcov?
              v <- vcov(object, ...)
              ses <- sqrt(diag(v))[parm]
              ci[] <- cf[parm] + ses %o% fac
              ci <- cbind(Lag = lags, ci, Estimate = cf[parm])
              if(se)
                  ci <- cbind(ci, StdError = ses)
              ci
          }
          )

.plot.acf.test <- function (x, y, data, method = "LiMcLeod", nlags, interval = 0.95,
                            legend = TRUE, ylim = c(NA, NA), ylim.fac = 1.2,
                            xlab = "Lag",
			    ci.lty = 2:3, # 2017-08-29 new, corresponding changes below
                            ...){
    if(missing(nlags))
        nlags <- seq(from = maxLag(x), to = 1, by = -5)

    x.iid <- whiteNoiseTest(x, h0 = "iid", n = x@n, interval = interval, # x = data,
                            method = method, nlags = nlags)
    ci <- x.iid$ci
    y <- cbind(1:nrow(ci), x.iid$ci, Estimate = x[1:nrow(ci)])

    ylim.loc <- ylim
    if(is.na(ylim[1])) ylim.loc[1] <- min(y[, 2])
    if(is.na(ylim[2])) ylim.loc[2] <- max(y[, 3])

    if(!missing(data)){
        x.garch <- whiteNoiseTest(x, h0 = "garch", x = data, nlags = nlags,
                                     interval = interval)
        ylim <- c(min(ylim[1], x.garch$ci[ , 1]),
                  max(ylim[2], x.garch$ci[ , 2])
                  )
        ## modify 'ylim.loc' if the corresponding 'ylim' element is NA
        if(is.na(ylim[1])) ylim.loc[1] <- min(ylim.loc[1], x.garch$ci[ , 1])
        if(is.na(ylim[2])) ylim.loc[2] <- max(ylim.loc[2], x.garch$ci[ , 2])
    }else{ # todo: modify some parameters for this case.
        legend <- FALSE
    }

    ## replace NA's in 'ylim' with the computed values
    if(is.na(ylim[1])) ylim[1] <- ylim.fac * ylim.loc[1]
    if(is.na(ylim[2])) ylim[2] <- ylim.fac * ylim.loc[2]

    .ciplot(x = y,
            ## lag_0 = TRUE,
            ylim = ylim,
            ci.col = "brown",
            ci.lty = ci.lty[1],
            xlab = xlab, ## xlab = colnames(x)[1],
            ## ylab = "Estimate & CI",
            ## ann = TRUE,
            ...
            )

    if(!missing(data)){
        x.garch <- whiteNoiseTest(x, h0 = "garch", x = data, nlags = nlags)
        .cilines(1:nrow(x.garch$ci), x.garch$ci, col = "blue", lty = ci.lty[2])
    }

    ## see ?legend, example thanks to Uwe Ligges
    if(legend){
        legend("topright", legend = c("H0: iid", "H0: garch"),
               col = c("brown ", "blue"), lty = ci.lty )
    }


    invisible()
}

.cilines <- function(x, y, type = "l", col = "blue", lty = 2, ...){
    lines(x, y[ , 1], type = type, col = col, lty = lty, ...)
    lines(x, y[ , 2], type = type, col = col, lty = lty, ...)
    invisible()
}

.plot_acr <- function(x, y, xlim = c(0,max(x)), ylim = c(-1, 1), type = "h", ...){

    plot(x, y, xlim = xlim, ylim = ylim, type = type, # xlab = xlab, ylab = ylab,
         ...)

    abline(h=0)

    ## if(ann)
    ##     title(xlab = xlab, ylab = ylab, ...)
    ## else
    ##     title(...)

    invisible()

}

.ciplot <- function(x, y, lag_0 = TRUE, ylim = c(-1, 1),
                    ci.col = "blue", ci.lty = 2,
                    xlab = colnames(x)[1],
                    ylab = "Estimate & rejection levels",
                    ann = TRUE,
                    ...){
    lags <- xval <- x[ , 1]
    yval <- x[ , "Estimate"]
    ## lbval <- x[ , 2]
    ## ubval <- x[ , 3]

    ## TODO: need more care here.
    if(lag_0){
        xval <- c(0, xval)
        yval <- c(1, yval)
    }

    ## plot(numeric(0), xlim = c(0,max(xval)), ylim = ylim, xlab = xlab, ylab = ylab,
    ##      ann = FALSE,
    ##      ...)
    ## abline(h=0)
    ## lines(c(0,0), c(0,1))

    .plot_acr(xval, yval, xlim = c(0,max(xval)), ylim = ylim, xlab = xlab, ylab = ylab,
         ann = ann, ...)

    ## lines(xval, yval, type = "h")
    .cilines(lags, x[ , 2:3], type = "l", col = ci.col, lty = ci.lty)

    ## if(ann)
    ##     title(xlab = xlab, ylab = ylab, ...)
    ## else
    ##     title(...)

    invisible()
}

setMethod("plot", c(x = "SampleAutocorrelations", y = "matrix"),
          function (x, y, main = "Acf test", ...){

              .ciplot(x = y,
                    lag_0 = FALSE, ## ylim = c(-1, 1),
                    main = main,
                    ## ci.col = "blue", ci.lty = 2,
                    ## xlab = colnames(x)[1],
                    ## ylab = "Estimate & CI",
                    ## ann = TRUE,
                    ...)
              invisible()
          }
          )

setMethod("plot", c(x = "SampleAutocorrelations", y = "missing"),
          function (x, y,
                    # data, method = "LiMcLeod", nlags, legend = TRUE,
                    main = "Acf test",
                    ...){
              .plot.acf.test(x = x, main = main, ...)
          }
          )

setMethod("plot", c(x = "SamplePartialAutocorrelations", y = "missing"),
          function (x, y,
                    # data, method = "LiMcLeod", nlags, legend = TRUE,
                    main = "Pacf test",
                    ...){
              .plot.acf.test(x = x, main = main, ...)
          }
          )

acfOfSquaredArmaModel <- function(model, maxlag){
    ar <- if(is.null(model$ar)) numeric(0) else model$ar
    arsq <- - coef(polynom(c(1, - ar))^2)[-1]
    masq <-   coef(polynom(c(1,   model$ma))^2)[-1]
    Sg <- model$sigma2 ^ 2

    autocovariances(list(ar = arsq, ma = masq, sigma2 = Sg), maxlag = maxlag)
}

## naive implementation
nvcovOfAcf <- function(model, maxlag){
    R <- autocovariances(model, maxlag = maxlag)
    r <- R / R[0]
    Rg <- acfOfSquaredArmaModel(model, maxlag = 2 * maxlag)

    res <- matrix(NA_real_, nrow = maxlag, ncol = maxlag)
    for(k in 1:maxlag){
        for(l in k:maxlag){
            res[k, l] <- Rg[l - k] + Rg[l + k] - 2 * Rg[k] * r[l] - 2 * Rg[l] * r[k] + 2 * Rg[0] * r[k] * r[l]
            if(k != l)
                res[l, k] <- res[k, l]
        }
    }
    res / R[0]^2
}

nvcovOfAcfBD <- function(acf, ma, maxlag){
    ## eq. 7.2.6. BD, p. 222
    res <- matrix(NA_real_, nrow = maxlag, ncol = maxlag)
    for(j in 1:maxlag){
        for(i in j:maxlag){
            maxk <- ma + max(i, ma)
            val <- 0
            for(k in 1 : maxk){
                val <- val + (acf[k + i] + acf[abs(k - i)] - 2 * acf[i] * acf[k]) *
                             (acf[k + j] + acf[abs(k - j)] - 2 * acf[j] * acf[k])
            }
            res[i, j] <- val
            if(i != j)
                res[j, i] <- res[i, j]
        }
    }
    res
}

rgarch1p1 <- function(n, alpha, beta, omega, n.skip = 100){
   Esigma2 <- omega / (1 - alpha - beta)

   ntot <- n.skip + n
   eta <- rnorm(ntot)
   sigma2 <- eps <- numeric(ntot)

   sigma2[1] <- Esigma2 # use E sigma_t^2 as initial value
   eps[1] <- sqrt(sigma2[1]) * eta[1]
   for(t in 2:ntot){
       sigma2[t] <- omega + alpha * eps[t-1]^2 + beta * sigma2[t-1]
       eps[t] <- sqrt(sigma2[t]) * eta[t]
   }

   eps[(n.skip + 1):ntot]
}

Try the sarima package in your browser

Any scripts or data that you put into this service are public.

sarima documentation built on Aug. 23, 2018, 9:03 a.m.