R/sarima.R

In GeoBosh/sarima: Simulation and Prediction with Seasonal ARIMA Models

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

coreXarmaFilter <- function(x, eps, ar = numeric(0), ma = numeric(0),
                            p = length(ar), q = length(ma), n = length(x),
                            from = max(p, q) + 1, intercept = 0){
    ceps <- if(q > 0){
                filter(eps, c(1, ma), sides = 1L)    # prepend '1' to the coef's.
                    # 2016-11-20 commenting out since:
                    #     (a) filter() sets ceps[1:q] to NA (but this means that
                    #         if p > 0 and p < q and from < q the result will be NA's,
                    #         which is likely to be a programming error).
                    #     (b) if from > q, ceps[1:q] is ignored below, so these values don't
                    #         matter in this case.
                    # ceps[1:q] <- 0
            }else{
                eps
            }
                    # intercept can be a constant or a vector of same length as ceps (and x)
    ceps <- ceps + intercept

    if(p > 0){
        ## this assumes that from > p
        init <- x[(from - 1):(from - p)] # filter() requires init in reverse order;

                                       # for "recursive", don't prepend '1' to 'ar'
        ceps[from:n] <- filter(ceps[from:n], ar, method = "recursive", init = init)
    }

    x[from:n] <- ceps[from:n]
    x
}

xarmaFilter <- function(model, x = NULL, eps = NULL, from = NULL, whiten = FALSE,
                        xcenter = NULL, xintercept = NULL){
                         # 2014-02-01 changing to "[[" to avoid partial matching for "p", etc.
                         #            which was a long standing cause of puzzling errors.
    ar        <- model[["ar",        exact = TRUE]] # model$ar
    ma        <- model[["ma",        exact = TRUE]] # model$ma
    center    <- model[["center",    exact = TRUE]] # model$center
    intercept <- model[["intercept", exact = TRUE]] # model$intercept
       # TODO: replace the above with:
       #    list2env(model, envir = environment()) # but need to change some is.null() below
       # (no, it may be dangerous in a user facing function.)

    p <- length(ar)
    q <- length(ma)

    if(is.null(from))                             # set a value for "from" if not supplied
        from <- 1 + max(p,q)
    else if(from <= max(p,q))
        stop("'from' must be greater than max(p,q)")

    ct <- if(is.null(intercept)) 0  else intercept
    if(!is.null(xintercept))
        ct <- xintercept + ct


    if(is.null(eps)){ # set initial values of eps to zero.  copy 'x' to keep the attributes
                      # of x, eg if x is a "ts" object, and ensure consistency of returned
                      # value in the two cases for 'whiten'
        eps <- x
        eps[] <- 0    # makes sense mostly when whiten = TRUE
    }else if(is.null(x)){
        x <- eps
        x[] <- 0
    }

    if(is.null(center) && is.null(xcenter)){
        flag.center <- FALSE
        y <- x
    }else{
        flag.center <- TRUE
        mu <- if(is.null(xcenter)) 0 else xcenter
        if(!is.null(center))
            mu <-  mu + center

        y <- x - mu
    }

    n <- length(x)
    stopifnot(n > 0)

    if(whiten){ # compute residuals (whiten the series x if x is arma with these params)
        eps <- coreXarmaFilter(x = eps, eps = y, ar = - ma, ma = - ar,
                               p = q, q = p, n = n, from = from, intercept = - ct)
        eps
    }else{      # compute a series from given residuals ("colour"/unwhiten eps)
        y <- coreXarmaFilter(x = y, eps = eps, ar = ar, ma = ma,
                             p = p, q = q, n = n, from = from, intercept = ct )

        x[from:n] <- if(flag.center) (y+mu)[from:n]  # add back the center
                     else               y[from:n]
        x
    }
}

# 2017-05-26: (somewhat) incompatible change - now uses prepareSimSarima()
sim_sarima <- function(model, n = NA, rand.gen = rnorm,
                   n.start = NA, x, eps, xcenter = NULL, xintercept = NULL, ...){
    ## decide if this patch should be permanent.
    if(is.list(model) && is.null(model$sigma2))
        model$sigma2 <- 1

    fu <-  prepareSimSarima(model = model, x = x, eps = eps, n = n,
                     n.start = n.start, xintercept = xintercept,
                     rand.gen = rand.gen)
    res <- fu()

    ## TODO: 2020-03-02 sim_sarima is documented to return an object from class "ts"
    ##       but doesn't!
    ## eg:
    ## frequency <- (environment(fu))$model$nseasons
    ## start <- if(is.null(frequency) || frequency == 1) 1 else frequency
    ## res <- ts(res, start = start, frequency = frequency)
    ##
    ## But it may be better to add arguments for the class, start, end, ...

    res
}

model2filter <- function(model){
    if(!is(model, "VirtualFilterModel"))
       model <- do.call("new", c(list("SarimaModel"), model))

    co <-  modelCoef(model, "ArmaModel")
    n_ur <- nUnitRoots(model)

    list(ar = co$ar, ma = co$ma,      # get 'ar' and 'ma' from 'co'
         sigma2    = model@sigma2,    # the rest from 'model'; could use sigmaSq(model), etc.
         center    = model@center,
         intercept = model@intercept,
         n_ur = n_ur  # 2018-07-28 new component
         )
}

## helper function - ensure that the argument is a list with components
##        before, init and main
.beforeInitMain <- function(x){
    res <- list(before = numeric(0), init = numeric(0), main = numeric(0))

    if(missing(x) || is.null(x))
        res
    else if(inherits(x, "list")){
                        #  x[c("before", "init", "main")] would return elements named <NA>
                        #  for absent components
        if(!is.null(x$before)) res$before <- x$before
        if(!is.null(x$init)) res$init <- x$init
        if(!is.null(x$main)) res$main <- x$main
        res
    }else{
        res$main <- x
        res
    }
}

prepareSimSarima <- local({
    ## these variables will be shadowed by assignments in the body of prepareSimSarima
    ## they are defined for somewhat more clear code; 'R CMD check' is also happier.
    x <- eps <- NULL
    ar <- ma <- numeric(0)
    p <- q <- 0
    from <- NULL # TODO: rethink!

           # the default for n.start is NA, so that in future automatic choice may be offered
    function(model, x = NULL, eps = NULL, n, n.start = NA, xintercept = NULL,
             rand.gen = rnorm){

        eps <- .beforeInitMain(eps)
        x <- .beforeInitMain(x)
        ct <- .beforeInitMain(xintercept)

        if(length(eps$before) != length(x$before)){
            if(length(eps$before) == 0)
                eps$before <- numeric(length(x$before))
            else if(length(x$before) == 0)
                x$before <- numeric(length(eps$before))
            else
                stop("Lengths of xbefore and innovbefore must be equal if both are present.")
        }

        if(length(eps$init) != length(x$init)){
            if(length(eps$init) == 0)
                eps$init <- numeric(length(x$init))
            else if(length(x$init) == 0)
                x$init <- numeric(length(eps$init))
            else
                stop("Lengths of xinit and innovinit must be equal if both are present.")
        }

        n.before <- length(x$before)

        if(is.na(n.start))
            n.start <- 0

        if(missing(n) || is.null(n))
            n <- length(eps$main) + length(eps$init) - n.start
        else{ # adjust x$main and eps$main
            eps$main <- c(eps$main,  
                          numeric(n.start + n - length(eps$init) - length(eps$main)) )
            x$main <- c(x$main,
                          numeric(n.start + n - length(x$init) - length(x$main)) )
        }

        from <- 1 + length(x$before) + length(x$init)
        n.drop <- n.before + n.start

        ## TODO: use new variables further below; for now save the old values in eps0,x0
        eps0 <- eps
        x0 <- x

        eps <- c(eps$before, eps$init, eps$main)
        x   <- c(x$before, x$init, numeric(n.start + n - length(x$init)))
            ## instead, could use:
            ##     x <- c(x$before, x$init, x$main)
            ## but it is preferable to state explicitly that any values in x$main
            ## are discarded.

        ## this is guaranteed by the above code but state the intention:
        stopifnot(length(eps) == length(x))

        ctt <- if(length(ct$main) == 0)
                   0
               else if(length(ct$main) == 1){  # constant xintercept
                   ct$main
               }else if(length(ct$main) > 1){  # non-constant xintercept
                   if(length(ct$before) == 0) # TODO: NA on place of 0?
                       ct$before <- numeric(length(eps0$before))
                   if(length(ct$init) == 0)
                       ct$init <- numeric(length(eps0$init))

                   c(ct$before, ct$init, ct$main)
               }

        stopifnot(length(x) == length(eps),
                  length(ctt) == 1  ||  length(ctt) == length(eps)
                  )

        filtmodel <- model2filter(model)

        ar <- filtmodel$ar
        ma <- filtmodel$ma
        mu <- filtmodel$center
        intercept <- filtmodel$intercept

        sigma <- sqrt(filtmodel$sigma2)

        p <- length(ar)
        q <- length(ma)

        if(intercept != 0)
            ctt <- ctt + intercept

        flag.center <- mu != 0
        n0 <- n
        rand.gen0 <- match.fun(rand.gen)

        ## 2018-07-28 krapka - for the case when no unit roots and no initial values
        ##            => generate initial values      
        flag.stat.init.values <- filtmodel$n_ur == 0  && from == 1 # 2018-07-28
        mpq <- max(p,q)  ## TODO: what  if mpq = 0?
        if(flag.stat.init.values){
            acv <- autocovariances(filtmodel[c("ar", "ma", "sigma2")], maxlag = mpq)
            partial_sds <- sqrt(.comboAcvf(acv, "psigma2"))[-1]
            partial_coefs <- acf2AR(acv[])
        }

       

        ## TODO: tryabva po-palna proverka na razmernostite!

        ## TODO: drugi varianti, kontrolirani s argument or otherwise;
        ##      e.g., use separate functions instead of conditionals
        f <- function(n = n0, rand.gen, ...){
            if(missing(rand.gen))
                rand.gen <- rand.gen0
            nx <- length(x)


            ## simulate (standardised) innovations and multiply by sigma
            ## 2022-03-24 was: eps.main <- sigma * rand.gen(n.start + n, ...)
            eps <- if(from > 1) # 2018-07-28 added if()
                       c(eps[1:(from-1)], sigma * rand.gen(length(eps) - (from - 1), ...))
                   else
                       sigma * rand.gen(n.start + n, ...)

            ## TODO: may need to modify x if n != n0

            if(flag.stat.init.values && mpq > 0){ # TODO: needs more thought
                x[1:mpq] <- partial_sds * eps[1:mpq]
                if(mpq > 1){
                    for(i in 2:mpq){
                        x[i] <- x[i] + sum(x[1:(i-1)] * partial_coefs[i-1, 1:(i-1)])
                    }
                }
                x[1:mpq] <- x[1:mpq] + mu # since mu will be subtracted below
            }
            

            y <- if(flag.center)
                     x - mu
                 else
                     x

            stopifnot(length(x) == length(eps), mpq < length(x)) # 2022-03-24
            
            y <- coreXarmaFilter(x = y, eps = eps, ar = ar, ma = ma, p = p, q = q, # n = n,
                                 # from = from, TODO: temporary commenting out (v. 0.4-5)
                                 intercept = ctt )

            x[from:nx] <-
                if(flag.center)
                    (y + mu)[from:nx]
                else
                    y[from:nx]

#browser()
            if(n.drop > 0) # 2018-07-28 added if()
                x[-(1:n.drop)] # TODO: keep the class of x?
            else
                x              # TODO: but this branch will keep it!
            
        }
        class(f) <- "simSarimaFun"
        f
    }
})

print.simSarimaFun <- function(x, ...){
    print("A function for simulation of the following Sarima model.")

    e <- environment(x)

    print(e$model) # TODO: do this more friendly.

    print("The parameters of the simulation are:")
    print("TODO: this print method is unfinished.")

    invisible(NULL)
}

sarima.f <- function(past = numeric(length(ar)),
                     n = max(2*length(past),12),
                     ar = numeric(0), ma = numeric(0),
                     intercept = 0,
                     pasteps = numeric(length(ma)), trend = numeric(n)){
    p <- length(ar)
    q <- length(ma)
    m <- max(p,q)

    if(length(past) > m) # keep only m past values
        past <- past[length(past) + (-m+1):0]

    res <- c(numeric(max(m-p,0)), past, numeric(n))

    if(length(pasteps) > m) # keep only m past values
        pasteps <- past[length(pasteps) + (-m+1):0]
    eps <- c(numeric(max(m-q,0)), pasteps, numeric(n))  # more care needed


    trend <- c(numeric(m), trend)

    for(i in  m + (1:n) ){
        res[i] <- intercept + trend[i] + eps[i]
        if(p>0)
            for(j in 1:p){
                res[i] <- res[i] + ar[j]*res[i-j]
            }
        if(q>0)
            for(j in 1:q){
                res[i] <- res[i] + ma[j]*eps[i-j]
            }
    }

    res[-(1:m)]             # Return the forecasts only.
}

fun.forecast <- function(  past
                         , n=max(2*length(past),12)
                         #, trend = numeric(n)
                         , eps = numeric(n)
                         , pasteps
                         , ...
                         ){
    ## model <- sarima.mod(...)
    sarima <- new("SarimaModel", ...)

    co <- filterCoef(sarima, "BD")
    ar = co$ar
    ma = co$ma
    center = sarima@center
    intercept = sarima@intercept

    if(missing(past))
        past <-  numeric(length(ar))

    if(missing(pasteps))
        pasteps <-  numeric(length(ma))

    fullintercept <- intercept
    if(center != 0 && sarima@iorder == 0 && sarima@siorder == 0) # Why only in this case?
                                       # (ans: because otherwise (1 - sum(ar) ) = 0)
        fullintercept <- intercept + (1-sum(ar)) * center

    res <- sarima.f(past = past, n = n, ar = ar, ma = ma
                  , intercept = fullintercept
                  , pasteps = pasteps
                  )                                                # trend=
    ts(res)         # Return the forecasts only.  TODO: set period?
}

## ## for now use sim_sarima as default
## simSarima <- sim_sarima
##
## setGeneric("simSarima", signature = c("model"))
##
## 2016-10-26 commenting out temporarily after removal of arCoef and maCoef
##
## setMethod("simSarima", "ArmaModel",
##           function(model, init = NULL, rand.gen = rnorm, info = "print", ...){
##               mo <- sarima.mod(ar = arCoef(model), ma = maCoef(model), mean = mean(model))
##               sim_sarima(model = mo, ...)
##           }
##           )