R/modelClasses.R

In sarima: Simulation and Prediction with Seasonal ARIMA Models

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

setGeneric("modelOrder",
           function(object, convention, ...){ standardGeneric("modelOrder") })
setGeneric("modelCoef" ,
           function(object, convention, component, ...){ standardGeneric("modelCoef") })
setGeneric("modelPoly",
           function(object, convention, ...){ standardGeneric("modelPoly") })
setGeneric("modelPolyCoef",
           function(object, convention, lag_0 = TRUE, ...){
               standardGeneric("modelPolyCoef")
           },
           signature = c("object", "convention")
           )

setGeneric("sigmaSq", def = function(object){ standardGeneric("sigmaSq") })
setGeneric("modelCenter", def = function(object){ standardGeneric("modelCenter") })
setGeneric("modelIntercept", def = function(object){ standardGeneric("modelIntercept") })

setGeneric("nUnitRoots", def = function(object){ standardGeneric("nUnitRoots") })
setGeneric("isStationaryModel",
           def = function(object){ standardGeneric("isStationaryModel") })

setClass("BJ", slots = c(dummy = "character"))
setClass("SP", slots = c(dummy = "character"))
setClass("BD", slots = c(dummy = "character"))

setClass("VirtualMeanModel",  contains = c("VIRTUAL"))

setClass("VirtualAutocovarianceModel",  contains = c("VIRTUAL"))
setClass("VirtualAutocorelationModel",
         contains = c("VirtualAutocovarianceModel", "VIRTUAL"))

setClass("VirtualPartialAutocovarianceModel", contains = "VIRTUAL")
setClass("VirtualPartialAutocorelationModel",
         contains = c("VirtualPartialAutocovarianceModel", "VIRTUAL"))

setClass("VirtualStationaryModel",
         contains = c("VirtualAutocovarianceModel", "VirtualMeanModel"))

setClass("VirtualWhiteNoiseModel", contains = c("VirtualStationaryModel", "VIRTUAL"))

setClass("VirtualFilterModel",  contains = c("VIRTUAL"))

  ## setMethod("sigmaSq", "VirtualFilterModel", function(object) object@sigma2)

setClass("VirtualArmaModel",
         contains = c("VirtualFilterModel", "VirtualStationaryModel", "VIRTUAL")
         ## , prototype = list(acv = new("Arma"), mean = numeric(0))
        )
setClass("VirtualArModel", contains = c("VirtualArmaModel", "VIRTUAL") )
setClass("VirtualMaModel", contains = c("VirtualArmaModel", "VIRTUAL") )

setClass("VirtualIntegratedModel", contains = c("VirtualFilterModel", "VIRTUAL") )

setClass("VirtualSarimaModel", contains = c("VirtualIntegratedModel", "VIRTUAL") )

setClass("VirtualArimaModel", contains = c("VirtualSarimaModel", "VIRTUAL") )
setClass("VirtualAriModel", contains   = c("VirtualArimaModel", "VIRTUAL") )
setClass("VirtualImaModel", contains   = c("VirtualArimaModel", "VIRTUAL") )

setClass("VirtualAutocovarianceSpec", contains = "VIRTUAL",
         slots = c(acvf = "ANY")
         )

setClass("AutocovarianceSpec",
         slots = c(acvf = "numeric")
         )

## setClass("AutocorrelationSpec", contains = "???")

setClass("InterceptSpec",
         slots = c(center =  "numeric", intercept = "numeric", sigma2 = "numeric"),
         prototype = list(center = 0, intercept = 0, sigma2 = NA_real_)
         )

setMethod("sigmaSq", "InterceptSpec", function(object) object@sigma2)
setMethod("modelCenter", "InterceptSpec", function(object) object@center)
setMethod("modelIntercept", "InterceptSpec", function(object) object@intercept)

setClass("ArmaSpec", contains = c("ArmaFilter", "InterceptSpec"))

## setMethod("innovationVariances", "ArmaSpec", function(object) object@sigma2 )
## TODO: Mean, Intercept - ne, tezi may sa za "ArmaModel"

setMethod("initialize", "ArmaSpec",
          function(.Object, ..., ar, ma, mean, check = TRUE){ # .Object <- callNextMethod()
              .Object <- callNextMethod(.Object, ...)
              if(!missing(mean)){
                  if(.Object@center == 0 && .Object@intercept == 0)
                      .Object@center <- mean
                  else
                      ## Could do more here but hardly worth the effort
                      ## TODO: Document this behaviour.
                      stop(paste0("Use argument 'mean' only when 'center' and 'intercept' ",
                                  "are missing or zero"))
              }

              if(!missing(ar)) .Object@ar <- as(ar, "BJFilter")
              if(!missing(ma)) .Object@ma <- as(ma, "SPFilter")

              ## 2018--05-30 Commenting out since there are no slots sar and sma here
              ## if(!missing(sar)) .Object@sar <- as(sar, "BJFilter")
              ## if(!missing(sma)) .Object@sma <- as(sma, "SPFilter")

              ## TODO: check that  orders are > 0 before calling polynomial functions
              if(check){
                  if(any(abs(solve(filterPoly(.Object@ar))) <= 1))
                      warning("The AR polynomial is not stable.")
                  
                  ## for now, do not warn for roots on the unit circle in MA models.
                  ## TODO: this needs further work.
                  if(any(abs(solve(filterPoly(.Object@ma))) < 1))
                      warning("The model is not invertible.")
              }

              .Object
          }
          )

setClass("AutocovarianceModel",
         contains = c("VirtualAutocovarianceModel", "AutocovarianceSpec")
        )

## TODO: are these needed?
##
##      "PartialAutocovarianceModel"
##
## setClass("AutocorrelationModel",
##          contains = c("VirtualAutocorrelationModel", "AutocorrelationSpec")
##         )
##
## setClass("PartialAutocorrelationModel",
##          contains = c("VirtualPartialAutocorrelationModel", "PartialAutocorrelationSpec")
##         )

setClass("ArmaModel",
         contains = c("ArmaSpec", "VirtualArmaModel")
        )

## TODO: Check if this convoluted inheritance is fine.
setClass("ArModel", contains = c("VirtualArModel", "ArmaModel"))
setClass("MaModel", contains = c("VirtualMaModel", "ArmaModel"))

setMethod("initialize", "ArModel",
          function(.Object, ...){
              .Object <- callNextMethod()
              validObject(.Object)
              .Object
          }
          )
setMethod("initialize", "MaModel",
          function(.Object, ...){
              .Object <- callNextMethod()
              validObject(.Object)
              .Object
          }
          )

ArmaModel <- function(...){
    new("ArmaModel", ...)
}

ArModel <- function(...){
    new("ArModel", ...)
}

MaModel <- function(...){
    new("MaModel", ...)
}

setAs("ArmaModel", "ArModel",
      function(from){
          if(modelOrder(from)$ma > 0)
              stop("Cannot convert 'model' to AR since it contains MA terms")
          ## the rest of the code is fom the default generated by package "methods"
          obj <- new("ArModel")      
          as(obj, "ArmaModel") <- from
          obj
      })

setAs("ArmaModel", "MaModel",
      function(from){
          if(modelOrder(from)$ar > 0)
              stop("Cannot convert 'model' to MA since it contains AR terms")
          ## the rest of the code is fom the default generated by package "methods"
          obj <- new("MaModel")      
          as(obj, "ArmaModel") <- from
          obj
      })

setValidity("ArModel", function(object){
    if(object@ma@order > 0)
        "Moving average terms found in ArModel object."
    else
        TRUE
})

setValidity("MaModel", function(object){
    if(object@ar@order > 0)
        "Autoregressive terms found in MaModel object."
    else
        TRUE
})

setClass("SarimaSpec", contains = c("SarimaFilter", "InterceptSpec"))

setClass("SarimaModel", contains = c("VirtualSarimaModel", "SarimaSpec"))

## TODO: subset seasonal arima (only some roots of one, specify them by
##           sequence number, i.e. as k in  exp(2pi k/s) , k= 0, ..., s-1

setMethod("nUnitRoots", "SarimaSpec",
          function(object){
              if(is.na(object@nseasons))
                  object@iorder
              else
                  object@iorder + object@nseasons * object@siorder
          }
          )
setMethod("isStationaryModel", "SarimaSpec",
          function(object){
              object@iorder == 0  && object@siorder == 0
          }
          )

setMethod("isStationaryModel", "VirtualIntegratedModel",
          function(object){
              nUnitRoots(object) == 0
          }
          )

setMethod("nUnitRoots", "VirtualStationaryModel", function(object){ 0 })
setMethod("isStationaryModel", "VirtualStationaryModel", function(object){ TRUE })

setMethod("modelPoly", c("VirtualMonicFilter", "missing"),
          function(object){
              filterPoly(object)
         }
         )

setMethod("modelPoly", c("VirtualFilterModel", "character"),
          function(object, convention){
              if(class(object) == convention)
                  modelPoly(object)
              else{
                  convention <- new(convention)  # classes that allow other things here
                                                 # should provide their own methods
                  modelPoly(object, convention  = convention )
              }
          }
          )

setMethod("modelPolyCoef", c("VirtualMonicFilter", "missing"),
          function(object, lag_0 = TRUE){
              filterPolyCoef(object, lag_0 = lag_0)
          }
          )

setMethod("modelPolyCoef", c("VirtualFilterModel", "character"),
          function(object, convention){
              if(class(object) == convention)
                  modelPolyCoef(object)
              else{
                  convention <- new(convention)  # classes that allow other things here
                                                 # should provide their own methods
                  modelPolyCoef(object, convention  = convention )
              }
          }
          )

setMethod("modelPoly", c("SarimaModel", "ArmaFilter"),
          function(object, convention){
              wrk <- filterPoly(object)

              list(ar = wrk$fullarpoly,
                   ma = wrk$fullmapoly )
          }
          )

setMethod("modelPolyCoef", c("SarimaModel", "ArmaFilter"),
          function(object, convention, lag_0 = TRUE){
              wrk <- modelPoly(object, convention)

              if(lag_0)
                  list(ar = coef(wrk$ar), ma = coef(wrk$ma) )
              else
                  list(ar = coef(wrk$ar)[-1], ma = coef(wrk$ma)[-1] )
          }
          )

setMethod("modelCoef", c("VirtualFilterModel", "missing", "missing"),
          function(object){
              filterCoef(object)
          }
          )

setMethod("modelCoef", c("VirtualFilterModel", "character", "missing"),
          function(object, convention){
              if(class(object) == convention)
                  modelCoef(object)
              else{
                  convention <- new(convention)  # classes that allow other things here
                                                 # should provide their own methods
                  modelCoef(object, convention  = convention )
              }
          }
          )

## convenience method for ARMA with BJ convention
setMethod("modelCoef", c("VirtualFilterModel", "BJ", "missing"),
          function(object, convention){
              if(class(object)[1] == "ArmaModel")
                 filt <- modelCoef(object)
              else{
                 filt <- modelCoef(object, convention  = "ArmaModel" )
              }
              list(ar = filt$ar, ma = -filt$ma)
          }
          )
## convenience method for ARMA with "SP" convention
setMethod("modelCoef", c("VirtualFilterModel", "SP", "missing"),
          function(object, convention){
              if(class(object)[1] == "ArmaModel")
                 filt <- modelCoef(object)
              else{
                 filt <- modelCoef(object, convention  = "ArmaModel" )
              }
              list(ar = - filt$ar, ma = filt$ma)
          }
          )
## convenience method for ARMA with "BD" convention
setMethod("modelCoef", c("VirtualFilterModel", "BD", "missing"),
          function(object, convention){
              if(class(object)[1] == "ArmaModel")
                 modelCoef(object)
              else{
                     # 2018-05-30 was:
                     #     modelCoef(object, convention  = "ArmaModel" )
                     # TODO: needs careful check - include tests if there are none yet.
                  modelCoef(object, convention  = "ArmaFilter" )
              }
          }
          )

setMethod("modelCoef", c("ArmaModel", "ArmaFilter", "missing"),
          function(object, convention){
              filterCoef(object)
      }
      )

setMethod("modelCoef", c("SarimaModel", "SarimaFilter", "missing"),
          function(object, convention){
              filterCoef(object)  # since "SarimaModel" inherits from "Sarima"Filter"
      }
      )

setMethod("modelCoef", c("SarimaModel", "ArmaFilter", "missing"),
          function(object, convention){
              wrk <- filterPolyCoef(object, lag_0 = FALSE)
              list(ar = - wrk$fullarpoly, ma = wrk$fullmapoly)
          }
          )

setMethod("modelCoef", c("SarimaModel", "ArFilter", "missing"),
          function(object, convention){
              wrk <- filterPolyCoef(object, lag_0 = FALSE)
              if(length(wrk$fullmapoly) > 0)
                   stop("Model not Ar-like (has non-trivial moving average part)")

                   # note: currently the result is ArmaFilter, not ArFilter
              list(ar = - wrk$fullarpoly, ma = numeric(0))
          }
          )

setMethod("modelCoef", c("SarimaModel", "MaFilter", "missing"),
          function(object, convention){
              wrk <- filterPolyCoef(object, lag_0 = FALSE)
              if(length(wrk$fullarpoly) > 0)
                  stop("Model not MA-like (has non-trivial autoregressive part)")

              list(ar = numeric(0), ma = wrk$fullmapoly)
          }
          )

setMethod("modelCoef", c("SarimaModel", "ArModel", "missing"),
          function(object, convention){
              wrk <- filterPolyCoef(object, lag_0 = FALSE)
              if(length(wrk$fullmapoly) > 0)
                   stop("Model not Ar-like (has non-trivial moving average part)")

              list(ar = - wrk$fullarpoly, ma = numeric(0))
          }
          )

setMethod("modelCoef", c("SarimaModel", "MaModel", "missing"),
          function(object, convention){
              wrk <- filterPolyCoef(object, lag_0 = FALSE)
              if(length(wrk$fullarpoly) > 0)
                  stop("Model not MA-like (has non-trivial autoregressive part)")

              ## 2020-02-22 was: list(ar = numeric(0), ma = coef(wrk$fullmapoly)[-1])
              list(ar = numeric(0), ma = wrk$fullmapoly)
          }
          )

setMethod("modelOrder", c("VirtualFilterModel", "missing"),
          function(object){
              filterOrder(object)
          }
          )


setMethod("modelOrder", c("VirtualFilterModel", "character"),
          function(object, convention){
              if(class(object) == convention)
                  modelOrder(object)
              else{
                  convention <- new(convention)  # classes that allow other things here
                                                 # should provide their own methods
                  modelOrder(object, convention  = convention )
              }
          }
          )

setMethod("modelOrder", c("ArmaModel", "ArFilter"),
          function(object, convention){
              wrk <- filterOrder(object)
              if(wrk$ma != 0)
                  stop("Non-zero moving average order")
              wrk
          }
          )

setMethod("modelOrder", c("ArmaModel", "MaFilter"),
          function(object, convention){
              wrk <- filterOrder(object)
              if(wrk$ar != 0)
                  stop("Non-zero autoregressive order")
              wrk
          }
          )

## check that these methods are inherited for convention = "ArmaModel", "ArModel" etc.
## But: even if they are, they may need to be redefined (different meaning!, see below)
setMethod("modelOrder", c("SarimaModel", "ArmaFilter"),
          function(object, convention){
               wrk <- modelOrder(object)
               if(is.na(wrk$nseasons))
                   wrk$nseasons <- 0

               with(wrk, list(ar = ar + iorder + (sar + siorder) * nseasons,
                              ma = ma + sma * nseasons ) )
          }
          )

setMethod("modelOrder", c("SarimaModel", "ArFilter"),
          function(object, convention){
               wrk <- modelOrder(object, "ArmaFilter")
               if(wrk$ma != 0)
                   stop("Non-zero moving average order")
               wrk
          }
          )

setMethod("modelOrder", c("SarimaModel", "MaFilter"),
          function(object, convention){
               wrk <- modelOrder(object, "ArmaFilter")
               if(wrk$ar != 0)
                   stop("Non-zero autoregressive order")
               wrk
          }
          )

## !!! notice the difference to c("SarimaModel", "ArmaFilter")
setMethod("modelOrder", c("SarimaModel", "ArmaModel"),
          function(object, convention){
               wrk <- modelOrder(object)
               if(is.na(wrk$nseasons))
                   wrk$nseasons <- 0

               with(wrk, {stopifnot(iorder == 0, siorder == 0)
                          list(ar = ar + sar * nseasons,
                               ma = ma + sma * nseasons )
                         })
           }
           )

setMethod("modelOrder", c("SarimaModel", "ArModel"),
          function(object, convention){
               wrk <- modelOrder(object, "ArmaModel")
               if(wrk$ma != 0)
                   stop("Non-zero moving average order")
               wrk
          }
          )

setMethod("modelOrder", c("SarimaModel", "MaModel"),
          function(object, convention){
               wrk <- modelOrder(object, "ArmaModel")
               if(wrk$ar != 0)
                   stop("Non-zero autoregressive order")
               wrk
          }
          )

.slots2list <- function(object){
    nams <- slotNames(object)
    res <- lapply(nams, function(x) slot(object, x))
    names(res) <- nams
    res
}

setAs("SarimaModel", "list",
      function(from){
          res <- .slots2list(from)

          res$ar <- filterCoef(res$ar)
          res$sar <- filterCoef(res$sar)
          res$ma <- filterCoef(res$ma)
          res$sma <- filterCoef(res$sma)

          res
      }
      )

as.list.SarimaModel <- function(x, ...){
    as(x, "list")
}

setAs("ArmaSpec", "list",
      function(from){
          res <- .slots2list(from)

          res$ar <- filterCoef(res$ar)
          res$ma <- filterCoef(res$ma)

          res
      }
      )

as.list.ArmaModel <- function(x, ...){ as(x, "list") }
as.list.ArmaSpec <- function(x, ...){ as(x, "list") }

.copy_cis <- function(from, to){
    to@center <- from@center
    to@intercept <- from@intercept
    to@sigma2 <- from@sigma2
    to
}

setAs("SarimaFilter", "ArmaFilter",
      function(from){
          ## was: modelCoef(from, "ArmaFilter") # TODOOOOOO! - DONE!
	  filt <- filterPolyCoef(from, lag_0 = FALSE)
          new("ArmaFilter", ar = filt$fullarpoly, ma = filt$fullmapoly)
      })

setAs("VirtualSarimaModel", "ArmaModel",
      function(from){
          if(!isStationaryModel(from))
              stop("This SARIMA model is not stationary.")
          filt <- as(from, "ArmaFilter") # multiply out the polynomials
          to <- new("ArmaModel", filt)
          to <- .copy_cis(from, to)
          to
      })

## utility for use in printing.
## cis stands for Center, Intercept and Sigma2.
.print_cis <- function(object, unconditional = FALSE){
    intercept <- modelIntercept(object)
    center <- modelCenter(object)

    if(unconditional  ||  center != 0)
        cat("Center: ", center, "\n")
    if(unconditional  ||  intercept != 0  || center == 0)
        cat("Intercept: ", intercept, "\n")

    cat("SigmaSq: ", object@sigma2, "\n")
}

## mis stands for Mean, Intercept and Sigma2.
.print_mis <- function(object){
    intercept <- modelIntercept(object)
    center <- modelCenter(object)

    if(is.na(intercept) || is.na(center))
        cat("mean: ", NA, "\n")
    else{##both not NA below
        pofzero <- sum(c(1, filterPolyCoef(object)$ar))
        fullintercept <- intercept +  pofzero * center # TODO: guard against div by 0
        mean <- fullintercept / pofzero
        if(intercept == 0)
            cat("mean: ", center, "\n")
        else if(center == 0)
            cat("intercept: ", intercept, "\n")
        else{ #both non-zero
            cat("mean: ", mean, "\n")
            cat("intercept: ", intercept, "(full intercept: ", fullintercept, ")", "\n")
        }
    }
    cat("sigmaSq: ", sigmaSq(object), "\n")
}

.print_formula <- function(object){
    order <- modelOrder(object)

    ar <- sar <- ma <- sma <- d <- ds <- intercept <- ""

    if(order$iorder > 0){
        d <- "(1-B)"
        if(order$iorder > 1)
            d <- paste0(d, "^", order$iorder)
    }

    if(order$siorder > 0){
        ds <- "(1-B^s)"
        if(order$siorder > 1)
            ds <- paste0(ds, "^", order$siorder)
    }

    if(order$ar > 0)  ar <- "Phi(B)"
    if(order$sar > 0) sar <- "Phi_s(B^s)"

    if(order$ma > 0)  ma <- "Theta(B)"
    if(order$sma > 0) sma <- "Theta_s(B^s)"

    arall <- paste0(d, ds, ar, sar)
    maall <- paste0(ma, sma)

    x <- "X(t)"
    ## TODO: take care of NA's
    if(object@center != 0){
        x <- paste0(x, " - ", "center")
        if(nchar(arall) > 0)
            x <- paste0("(", x, ")")
    }

    intercept <- if(object@intercept != 0)
                    "intercept + "

    e <- "e(t)"
    cat("Model: ", arall, x, " = ", intercept, maall, e, "\n", sep = "")
}

setMethod("show",
    signature(object = "InterceptSpec"),
    function (object){
        .reportClassName(object, "InterceptSpec")
        .print_cis(object, unconditional = TRUE)
        invisible(NULL)
    }
)

## TODO: methods for these!
##
## if(model$mean != 0 && d==0 && ds==0)
##    res$fullintercept <- res$intercept + (1-sum(res$fullar))*model$mean
##
##  res$fullmean <- res$mean                                # more care needed!
##  if(model$intercept != 0 && d==0 && ds==0)
##    res$fullmean <- res$fullintercept/(1-sum(res$fullar))

setMethod("show",
    signature(object = "ArmaModel"),
    function (object){
        .reportClassName(object, "ArmaModel")
        .print_mis(object)
        callNextMethod()
        invisible(NULL)
    }
)
setMethod("show",
    signature(object = "ArModel"),
    function (object){
        .reportClassName(object, "ArModel")
        ## TODO: avoid printing the ma slot?
        callNextMethod()
    }
)
setMethod("show",
    signature(object = "MaModel"),
    function (object)
    {
        .reportClassName(object, "MaModel")
        ## TODO: avoid printing the ar slot?
        callNextMethod()
    }
)

setMethod("show",
    signature(object = "SarimaModel"),
    function (object){
        .reportClassName(object, "SarimaModel")
        .print_formula(object)
        cat("\n")
        if(isStationaryModel(object))
            .print_mis(object)
        else
            .print_cis(object)
        callNextMethod()
    }
)

summary.SarimaSpec <- function(object, ...){
    cat("Intercept: ", object@intercept, "\n")
    cat("Center: ", object@center, "\n")
    cat("Innovation variance: ", object@sigma2, "\n")

    summary.SarimaFilter(object)
}

summary.SarimaModel <- function(object, ...){
    .reportClassName(object, "SarimaModel")
    summary.SarimaSpec(object)
}

as.SarimaModel <- function (x, ...) UseMethod("as.SarimaModel")

as.SarimaModel.Arima <- function (x, ...){
    order <- x$arma
    p <- order[1]
    q <- order[2]
    P <- order[3]
    Q <- order[4]
    nseasons <- order[5]
    d <- order[6]
    D <- order[7]

    ## TODO: include the mean (named 'intercept' in Arima objects)
    co <- coef(x)
    ar <- co[seq_len(p)]
    ma <- co[p + seq_len(q)]
    sar <- co[p + q + seq_len(P)]
    sma <- co[p + q + P + seq_len(Q)]
    
    new("SarimaModel", ar = ar, ma = ma, sigma2 = x$sigma2,
                       sar = sar, sma = sma, nseasons = nseasons,
                       iorder = d, siorder = D)
}

## setClass("Spectrum", slots = c(freq = "numeric", spec = "numeric", model = "ANY"),
##          contains = "function")

setClass("Spectrum", slots = c(call = "call", model = "ANY"), contains = "function")
setClass("ArmaSpectrum", contains = "Spectrum")

setMethod("initialize", "ArmaSpectrum",
          function(.Object, ar = numeric(0), ma = numeric(0), sigma2 = NA_real_, ...){
              .Object <- callNextMethod(.Object, ...)
              if(is.null(.Object@model))
                  .Object@model <- ArmaModel(ar = ar, ma = ma, sigma2 = sigma2)
              ## else
              ##   something wrong here, easy to get inconsistencies,
              ##   Don't document for users!
              
              if(is.null(body(.Object@.Data))){ # function not set, do it
                  e <- new.env(parent = topenv())
                  e$ar <- ar
                  e$ma <- ma
                  e$sigma2 <- sigma2
                  freq <- standardize <- NULL ## otherwise 'R CMD check' complains
                  f <- as.function(alist(freq = , standardize = TRUE,
                                         .spdArma(ar = ar, ma = ma, sigma2 = sigma2,
                                                  freq = freq, standardize = standardize)$spec
                                         ), envir = e)
                  .Object@.Data <- f
              }
              .Object
          })

plot.Spectrum <- function(x, y, to, from = y, n = 128, standardize = TRUE,
                          log = NULL, main = "Spectral density", xlab = "Frequency", 
                          ylab = NULL, ...){
    if(missing(from))
        from <- 0
    if(missing(to))
        to <- 0.5
    if(is.null(ylab)){
        ylab <- if(identical(log, "y"))
                    "log(Spd)"
                else
                    "Spd"
    }
    ## else leave as is

    if(standardize)
        plot.function(x, type = "l", from = from, to = to, n = n,
                      log = log, main = main, xlab = xlab, ylab = ylab, ...)
    else{
        ## f <- function(z) x(z, standardize = FALSE)
        ## curve(f, from = from, to = to, ...)
        curve(x(x, standardize = FALSE), from = from, to = to, n = n,
              log = log, main = main, xlab = xlab, ylab = ylab, ...)
    }
}

setMethod("plot", "Spectrum", plot.Spectrum)

format.Spectrum <- function (x, ..., n = 128, standardize = TRUE){
    freq <- seq(0, 0.5, length.out = n)
    spec <- x(freq)

    res <- ""

    res <- c(res, "\nPeaks:")
    max_flags <- .local_maxima(spec)
    ma_freq <- freq[max_flags]
    ma_spec <- spec[max_flags]
    pma <- 1 / ma_freq
    pma <- ifelse(is.infinite(pma), 0, pma)
    res <- c(res, capture.output(print(cbind(freq = ma_freq, spec = ma_spec, period = pma), ...)),
                  "")
    
    res <- c(res, "Troughs:\n")
    min_flags <- .local_minima(spec)
    mi_freq <- freq[min_flags]
    mi_spec <- spec[min_flags]
    pmi <- 1 / mi_freq
    pmi <- ifelse(is.infinite(pmi), 0, pmi)
    res <- c(res, capture.output(print(cbind(freq = mi_freq, spec = mi_spec, period = pmi), ...)),
                  "")

    p <- max(ma_spec) / min(mi_spec)
    res <- c(res, "max peak/min trough:")
    res <- c(res, paste0("\t", format(p, ...)))
    res <- c(res, "")

    res
}

format.ArmaSpectrum <- function (x, ..., n = 128, standardize = TRUE){
    ar     <- environment(x@.Data)$ar
    ma     <- environment(x@.Data)$ma
    sigma2 <- environment(x@.Data)$sigma2
    
    res <- c(paste0(if(standardize) "standardized " else "",
                    "spectral density of the following ",
                    "ARMA(", length(ar), ",", length(ma), ") model:"),
             paste0("  ar coef: ", 
                    paste0(formatC(ar, zero.print = TRUE, ...), collapse = " ")),
             paste0("  ma coef: ", 
                    paste0(formatC(ma, zero.print = TRUE, ...), collapse = " ")),
             paste0("  sigma2:  ", format(sigma2, ...)),
             "")

    if(standardize)
        sigma2 <- 1
    else if(is.na(sigma2))
        stop("sigma2 must be non-NA when 'standardize' is FALSE")

    if(length(ar) == 0 && length(ma) == 0){
        res <- c(res, paste0("constant, equal to ", sigma2, " for all frequencies"))
	return(res)
    }

    res <- c(res, format.Spectrum(x, ..., n = n, standardize = standardize))

    res
}

print.Spectrum <- function (x, ..., n = 128, standardize = TRUE){
    wrk <- format(x, ..., n = n, standardize = standardize)
    cat(wrk, sep = "\n")
    if(length(list(...)) == 0 && missing(n) && missing(standardize) )
        plot(x)

    invisible(x)
}

setMethod("show", "Spectrum", function(object){
    cat(format(object), sep = "\n")
    plot(object)
    invisible(NULL)
})