R/autocovariances.R

In sarima: Simulation and Prediction with Seasonal ARIMA Models

## 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")){
        switch(slot,
               varnames = {
                   if(length(slot(x, slot)) == 0){
                       cat('Slot ', slot, ':  ', " <not set>", "\n", sep = "")
                   }else{
                       cat('Slot ', slot, ':', "\n", sep = "")
                       print(slot(x, slot))
                   }
               },
               objectname = { 
                   cat('Slot ', slot, ':  ', slot(x, slot), "\n", sep = "")
               },
               ## default       
               { cat('Slot ', slot, ':', "\n", sep = "")
                   print(slot(x, slot))
               }
        )
    }
}

setMethod("show", signature(object = "Autocovariances"),
          function (object)
          {
              .reportClassName(object, "Autocovariances")
                  # 2019-03-29 (similarly below) was: callNextMethod()
                  #    is it better to drop prefix "Lag_":
                  #        dataWithLagNames(object, "") 
              print(dataWithLagNames(object)) 
          })

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

setMethod("show", signature(object = "PartialAutocovariances"),
          function (object)
          {
              .reportClassName(object, "PartialAutocovariances")
              print(dataWithLagNames(object)) 
          }
          )

setMethod("show", signature(object = "PartialAutocorrelations"),
          function (object)
          {
              .reportClassName(object, "PartialAutocorrelations")
              print(dataWithLagNames(object)) 
          }
          )

setMethod("show", signature(object = "PartialVariances"),
          function (object)
          {
              .reportClassName(object, "PartialVariances")
              print(dataWithLagNames(object)) 
          }
          )

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 = "SamplePartialAutocovariances"),
          function (object)
          {
              .reportClassName(object, "SamplePartialAutocovariances")
              callNextMethod()
              .printFittedPart(object)
          }
          )

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

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

.basecfclass2S4 <- function(x, r0){
    stopifnot(inherits(x, "acf")) # 2020-02-29 was: stopifnot(class(x) == "acf"), similarly below

    lagged <- Lagged(x) # Lagged knows how to interpret "acf" objects
    if(!missing(r0))
        lagged[0] <- if(inherits(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(res@data) == 3)
            #     new("FlexibleLagged", data = obj@data[ -1, , ])
            # else
            #     new("FlexibleLagged", data = obj@data[-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, ...)
            ## BUGFIX: was: res <- as(obj, "FlexibleLagged")
            res <- new("FlexibleLagged", data = pacr)
            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 <- if(missing(maxlag))
                   autocovariances(x, ...)
               else
                   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")

setMethod("autocovariances",
          signature(x = "VirtualAutocovariances"),
          function (x, ...)
          {
              stop("Can't compute autocovariances from objects from class ", class(x))
          }
          )

setMethod("autocovariances",
          signature(x = "Autocovariances", maxlag = "missing"),
          function (x, ...)
          {
              x
          }
          )

setMethod("autocovariances",
          signature(x = "Autocovariances", maxlag = "ANY"),
          function (x, maxlag, lag_0, ...)
          {
              x[] <- x[0:maxlag]   # may introduce NA's
              x
          }
          )

## 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(obj@data)) && dim(obj@data) == 3)
                #     obj@data[ -1, , ]
	        # else
                #     obj@data[-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 object from class ", class(x))

              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.')
              ## }
          }
          )

setMethod("autocorrelations",
          signature(x = "Autocorrelations", maxlag = "missing", lag_0 = "missing"),
          function (x, ...)
          {
              ## as(x, "Autocorrelations")   # x may be from a subclass
              x
          }
          )

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
              x[] <- x[0:maxlag]   # may introduce NA's
              x
          }
          )



setMethod("autocorrelations",
          signature(x = "Autocovariances", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              if(missing(maxlag))
                  maxlag <- maxLag(x)
              acr <- x[0:maxlag] / x[[0]] # TODO: Only univariate case for now!
              if(is(x, "Fitted"))
                  new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
              else
                  new("Autocorrelations", data = acr)
          }
          )

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

setMethod("autocorrelations",
          signature(x = "PartialAutocorrelations", maxlag = "ANY", lag_0 = "missing"),
          function (x, maxlag, ...)
          {
              acr <- modelCoef(x, "Autocorrelations")
              res <- if(is(x, "Fitted"))
                         new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
                     else
                         new("Autocorrelations", data = acr)
              if(!missing(maxlag))
                  ## BugFix 2020-02-29 was: res@data <- res[0:maxlag]
                  ##     TODO: needs more testing
                  res@data <- Lagged(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-stationary
        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)
              ## 2020-02-29 not needed, maxlag has already been used above!
              ##     Also, this is wrong since res@data has a Lagged class.
              ## if(!missing(maxlag))
              ##     res@data <- res[0:maxlag]     # may introduce NA's
              res
          }
          )


## BUGFIX: this cannot be done, see details elsewhere in this document!
## 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(acvf@data, 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.")
##       })

setAs("Autocorrelations", "Autocovariances",
      function(from){
          new("Autocovariances", data = from[])
      })

setAs("SampleAutocorrelations", "SampleAutocovariances",
      function(from){
          new("SampleAutocovariances", data = from[], as(from, "Fitted"))
      })

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("VirtualAutocovariances", "Autocovariances", "missing"),
          function(object, convention){
              stop("Can't obtain autocovariances from object from class ", class(object))
          }
          )

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){
              ## 2020-02-29 changing this (here and at similar places) with temporary patch
              ##     since Lagged2d objects don't have method for "character" (also "logical")
              ## 
              ## TODO: UPDATE PACKAGE "lagged" with such methods, if sensible!
              ##
              ##      object@data["acvf", ]
              tmp <- which(rownames(object@data[]) == "acvf")
              object@data[tmp, ]
          }
          )

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

setMethod("modelCoef", c("ComboAutocovariances", "PartialAutocovariances", "missing"),
          function(object, convention){
              ## object@data["pacvf", ]
              tmp <- which(rownames(object@data[]) == "pacvf")
              object@data[tmp, ]
          }
          )

setMethod("modelCoef", c("ComboAutocovariances", "PartialVariances", "missing"),
          function(object, convention){
              ## object@data["sigma2", ]
              tmp <- which(rownames(object@data[]) == "sigma2")
              object@data[tmp, ]
          }
          )


## 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", ]
              tmp <- which(rownames(object@data[]) == "acf")
              object@data[tmp, ]
          }
          )

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

setMethod("modelCoef", c("ComboAutocorrelations", "PartialAutocorrelations", "missing"),
          function(object, convention){
              ## object@data["pacf", ]
              tmp <- which(rownames(object@data[]) == "pacf")
              object@data[tmp, ]
          }
          )

## need class StandardizedPartialVariances
## setMethod("modelCoef", c("ComboAutocorrelations", "StandardizedPartialVariances", "missing"),
##           function(object, convention){
##           object@data["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 <- object@data
                  #     pacfonly <- pacf[-1]
              pacfonly <- object[1:maxLag(object)]
                   # 2022-02-14 was:  ar <- FitAR::PacfToAR(pacfonly)
              ar <- pacf2Ar(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
##               stop("Partial autocorrelations are not uniquely determined by partial autocovariances")
##           }
##           )
## 
## setMethod("modelCoef", c("PartialAutocovariances", "Autocorrelations", "missing"),
##           function(object, convention){
##               ## pacfonly <- object[] / object[0] # TODO: currently scalar only
##               ## ar <- FitAR::PacfToAR(pacfonly[-1]) # drop lag 0 !
##               ##      # alternatively:
##               ##      # acvf <- ltsa::tacvfARMA(phi = ar, maxLag = length(ar))
##               ##      # acf <- acvf/acvf[0]
##               ## ARMAacf(ar, lag.max = length(ar))
##               stop("Autocorrelations are not uniquely determined by partial autocovariances")
##           }
##           )

## 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, npar = npar, 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, npar = npar, nlags = nlags, method = method, interval = interval)
          })

acfMaTest <- function(acf, ma,  n, nlags,
                       interval = 0.95){
    if(!is(acf, "Lagged"))
        acf <- Lagged(acf)
    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)

        ## 2019-03-15 <BUGFIX> was:  int <- cbind(int, - int)
        ##         also need levels for the lags!
        rho_hat_null <- c(acf[seq_len(ma)], rep(0, length(int) - ma))
        int <- rho_hat_null + 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
}

acfWnTest <- function(acr, x, nlags, interval = 0.95, ...){
    maxlag <- max(nlags)
    n <- length(x)

    v <- nvarOfAcfKP(x, maxlag = maxlag, ...)
    ## nvcov is diagonal in this case
    Q <-  n * cumsum(acr[1:maxlag]^2 / v )[nlags]
    
    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))
        se <- sqrt(v / n)
        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
           },
           "sv" = ,
           "arch-type" = {
               acfWnTest(object, ...)
           },
           ##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
              if(inherits(assuming, "Arima")){
                  ## todo: could get vcov of params and transform it for acf's
                  ##       maybe have additional argument here (method?)
                  mo <- list(ar     = assuming$model$phi,
                             ma     = assuming$model$theta,
                             sigma2 = assuming$sigma2 )
                  res <- nvcovOfAcf(mo, maxlag = maxlag) / object@n
                  res
              }else
                  vcovAcf(object, assuming, maxlag = maxlag, ...)
          }
          )

setGeneric("vcovAcf", function(object, assuming, ...) stop("no suitable method"))

setMethod("vcovAcf", c("SampleAutocorrelations", "character"),
          function (object, assuming, maxlag = maxLag(object), ...) {
              if(length(assuming) > 1)
                  stop("'assuming' must be of length 1")
              
              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("vcovAcf", c("SampleAutocorrelations", "ArmaModel"),
          function (object, assuming, maxlag = maxLag(object), ...) {
              res <- nvcovOfAcf(c(modelCoef(assuming), list(sigma2 = sigmaSq(assuming))), 
                                maxlag = maxlag) / object@n
              res
          }
          )

se <- function (object, ...) {
    vc <- vcov(object, ...)
    sqrt(diag(vc))
}
## TODO: make S4 generic

.fixed_values_under_h0 <- function(object, lags){
    m <- length(lags)
    if(identical(object, "iid"))
        numeric(m)
    else if(identical(object, "garch"))
        numeric(m)
    else if(inherits(object, "Arima")){
        order <- object$arma
        names(order) <- c("ar", "ma", "sar", "sma", "period", "d", "ds")
        if(order["d"] != 0 || order["ds"] != 0)
            stop("Can't compute acf confint for non-stationary models")
        else if(all(order[c("ar", "sar")] == 0)){
            if(order["sma"] == 0){
                ## TODO: can be refined
                res <- ifelse(lags <= order["ma"], NA_real_, 0)
                res
            }else if(order["ma"] == 0){ # pure seasonal MA
                res <- ifelse(lags %% order["period"] == 0  &
                              lags %/% order["period"] <= order["ma"] ,
                              NA_real_, 0)
                res
            }else
                res <- rep(NA_real_, m)
        }
    }else if(is(object, "ArmaModel")){
        ## for a theoretical model all acfs are fixed
        ## for fitted models => see the case 'Arima'
        ## TODO: need further methods here
            # if(modelOrder(object)$ar == 0)
            #     ifelse(lags <= modelOrder(object)$ma, NA_real_, 0)
            # else
#browser()
        autocovariances(c(modelCoef(object), list(sigma2 = sigmaSq(object))), 
                        maxlag = max(lags))[lags]
    }else ## default method - all not fixed
        rep(NA_real_, m)
}

# setMethod("diagOfVcov", "SampleAutocorrelations",

setMethod("confint", "SampleAutocorrelations",
          function (object, parm, level = 0.95, se = FALSE, maxlag, ..., assuming){
              ## 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?
                  # 2022-01-17 was:
                     # v <- vcov(object, ...)
                     # ses <- sqrt(diag(v))[parm]
              ses <- se(object, assuming = assuming, ...)[lags]

              cg <- cf[parm]
              if(!is.null(assuming)){
                  int_types <- .fixed_values_under_h0(assuming, lags)
                  if(!all(is.na(int_types))){
                      cg[!is.na(int_types)] <- int_types[!is.na(int_types)]
                  }
              }else
                  int_types <- rep(NA_real_, length(parm))

              ci[] <- cg + ses %o% fac
              ci <- cbind(Lag = lags, ci, Estimate = cf[parm], H0_fixed = int_types)
              if(se)
                  ci <- cbind(ci, StdError = ses)
              attr(ci, "assuming") <- assuming
              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
                            h0 = "garch", # 2019-04-09 new to accommodate the new white noise test
                                          #                (Kokoszka-Politis)
                            ...){
    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 = h0, 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", paste0("H0: ", h0)),
               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 = NULL, maxlag){
    maacf <- new("Lagged1d", data = acf[])
    max_avail_lag <- maxLag(maacf)

    if(is.null(ma))
        ma <- max_avail_lag
    else
        stopifnot(ma <= max_avail_lag)

    ## TODO: this may make max_avail_lag too large. Fix!
    if(max_avail_lag < 2 * maxlag + ma)
        max_avail_lag <- 2 * maxlag + ma

    if(ma < max_avail_lag)
        maacf[(ma + 1) : max_avail_lag] <- 0

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

    res
}

nvarOfAcfKP <- function(x, maxlag, center = FALSE, acfscale = c("one", "mom")){
    acfscale <- match.arg(acfscale)
    if(center)
        x <- x - mean(x)
    acv <- acf(x^2, lag.max = maxlag, demean = FALSE, type = "covariance", plot = FALSE)
    acv <- acv$acf[ , , ] # make numeric, or more explicitly: acv$acf[ , , 1]
    ## the denominator, NOTE: it is the square of the lag 0 acv of x, not x^2
    n <- length(x)
    den <- (sum(x^2) / n) ^ 2
    switch(acfscale,
           one = acv[-1] / den,
           mom = {
               fac <- n / (n - seq_len(maxlag))
               fac * acv[-1] / den
           }
           ## 2020-02-29 removing this since match.arg() above catches this error.
           ## ,
           ## ## default
           ## stop("argument 'acfscale' must be one of 'one', 'mom' or abbreviation thereof")
           )
}

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]
}