R/fitTools.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.

## todo: expand_ar() seems unused, maybe it is obsolete(d)?
expand_ar <- function(ar, sar, s){
    x <- polynom(c(0,1))
    phi <- polynom(c(1, -ar))
    sphi <- polynom(c(1, -sar))(x^s)

    coef(phi * sphi)  
}

poly_x     <- polynom(c(0,1))
poly_delta <- polynom(c(1, -1)) # = 1 - poly_x

.sarima_env <- function(){

    .process_fixed <- function(x, fixed = NULL){
        if(is.null(x))
            stop("'x' must have a value here")
        xlen <- length(x)
        
        ## (2018-08-23) TODO: this if/else sequence needs cleaning up!
        if(is.null(fixed))
            rep(FALSE, xlen)
        else if(isTRUE(fixed))
            rep(TRUE, xlen)
        else if(is.logical(fixed)){
            ## 2018-08-23 fiixng missing case length=1 and FALSE
            ##   TODO: check nevertheless if this was omitted on purpose?
            if(length(fixed) == 1)
                rep(fixed, xlen)
            else if(length(fixed) == xlen)
                fixed
            else
                stop("if 'fixed' is logical, its length must be equal to one or the order.")
        }else if(is.numeric(fixed)){
            ## indices of fixed elements, must be whole numbers
            if(!is.integer(fixed) && any(round(fixed) != fixed))
                stop("if numeric, 'fixed' must contain integer numbers")
            if(any(fixed <= 0) || any(fixed > xlen))
                stop("if numeric, 'fixed' must contain integers between 1 and 'order'")
            res <- rep(FALSE, xlen)
            res[fixed] <- TRUE
            res # needed since 'res[fixed] <- value'  returns the rhs
        }else{
            stop("unsupported type of 'fixed'")
        }
    }

    ## currently essentially identical to .process_fixed()
    .process_atanh.tr <- function(x, transf = NULL){
        if(is.null(x))
            stop("'x' must have a value here")
        xlen <- length(x)
        
        if(is.null(transf))
            rep(FALSE, xlen)
        else if(isTRUE(transf))
            rep(TRUE, xlen)
        else if(is.logical(transf)){
            if(length(transf) == 1)
                rep(transf, xlen)
            else if(length(transf) == xlen)
                transf
            else
                stop("if 'transf' is logical, its length must be equal to one or the order.")
        }else if(is.numeric(transf)){
            ## indices of transf elements, must be whole numbers
            if(!is.integer(transf) && any(round(transf) != transf))
                stop("if numeric, 'transf' must contain integer numbers")
            if(any(transf <= 0) || any(transf > xlen))
                stop("if numeric, 'transf' must contain integers between 1 and 'order'")
            res <- rep(FALSE, xlen)
            res[transf] <- TRUE
            res
        }else{
            stop("unsupported type of 'transf'")
        }
    }

    ## parcor
    ## roots ?
    ## fixed
    ## nonfixed
    ## convention ?
    ## dispname ?
    ##

    ## based on co_ar() in sarima() and ar() in .sarima_env
    .set_coef <- function(order, coef = NULL, fixed = NULL, nonfixed = NULL,
                          atanh.tr = NULL, ...){
        wrk <- NULL
        coef.type <- "coef"  # todo: any subtleties?
        if(!is.null(order))
            wrk <- rep(NA_real_, order)
    
        if(!missing(coef) && !is.null(coef)){
            if(is.null(wrk))
                wrk <- coef
            else{
                    # todo: is it reasonable to allow a smaller length, so that coef's from
                    #       lower order fits can be used as initial values more easily?
                ## Note (JH, 11/10/18): I like that length(coef) == length(wrk)
                ##    but understand it is not necessary.
                ##    If this is true, 'p' is not required (see note below)
                if(length(coef) != length(wrk))
                    stop("if given, 'coef' must have length equal to 'order'")
                ## 2018-08-23 rewriting using seq_along:
                wrk[seq_along(coef)] <- coef  # wrk[1:length(coef)] <- coef
            }
        }
    
        a <- .process_fixed(wrk, fixed)
        b <- .process_atanh.tr(wrk, atanh.tr)

        if(!missing(nonfixed) || !is.null(nonfixed)){
                                # declare non-fixed some coef's previously marked as fixed
            ## !!! 2018-08-23 - this 'b' destroys the 'b' from atanh.tr, seems oversight!
            ##        Changing this 'b' to 'bb'.
            ##        TODO: verify nevertheless!
            ##   b <- .process_fixed(wrk, nonfixed)
            bb <- .process_fixed(wrk, nonfixed)
                # !bb is FALSE only for the elements declared nonfixed by 'nonfixed', so this
                # ! will transfer this to 'a', leaving others as they are (fixed or nonfixed)
            a <- a & !bb
        }

        res <- wrk
        attr(res, "fixed") <- a
        attr(res, "atanh.tr") <- b

        ## store addtional info, e.g. parcor, negate
        dots <- list(...)
        for(at in names(dots))
            attr(res, at) <- dots[[at]]
        
        res
    }

    ar <- function(p = 0, ar, sign = "-", atanh.tr = TRUE, ...){
        ## NOTE: do not name p and ar in the call for flexibility:
        coef <- .set_coef(p, ar, sign = sign, dispname = "ar", atanh.tr = atanh.tr, ...)
        list(name = "ar", p = p, coef = coef)    # TODO: remove 'p' from the return value?
        ## Note (JH, 11/10/18): 'p' is never used, I would agree on removing 'p'.
        ##   If one needed, p could be obtained from the length of coef
    }

    ma <- function(q = 0, ma, sign = "+", atanh.tr = TRUE, ...){
        coef <- .set_coef(q, ma, sign = sign, dispname = "ma", atanh.tr = atanh.tr, ...)
        list(name = "ma", q = q, coef = coef)
    }

    s <- function(...){
        ## TODO: for now 'm' of length 1
        ## Note (JH, 11/10/18): what else would be an argument here?
        dots <- c(...)
        if(length(dots) != 1)
            stop("currently only one argument is supported")
        p <- dots[1] - 1
        coef <- .set_coef(p, rep(1, p), sign = "+", fixed = TRUE, operator = TRUE, dispname = "s")
        list(name = "s", p = p, coef = coef)
    }

    ## unit roots
    i <- function(d, ...){
        ## TODO: this needs special method later!
        coef <- .set_coef(1, 1, sign = "-", fixed = TRUE, d = d, operator = TRUE,
                          dispname = "i", ...)
        list(name = "i", d = d, p = 1, coef = coef)
    }

    u <- function(u, fixed = TRUE, operator = all(fixed), ...){
                     # u may be a list to accommodate mixed complex and real values, probably
                     # not a good idea. Do not convert to use Recall() if u is 'list'
        f <- function(x){
            if(is.complex(x))
                -2*cos(Arg(x))    # equivalently: -2*Re(x), but this assumes that Mod(x)=1
            else -2 * cospi(2*x) # x is a fraction of 2pi
        }
        co <- sapply(u, f) 
        coef <- lapply(co, function(x) .set_coef(2, c(x, 1), sign = "+", fixed = fixed,
                                                 operator = operator, dispname = "su") )
        list(name = "u", u = u, coef = coef)
    }

    ## spec for unit roots to be estimated (unless fixed)
    ## TODO: require p >= 2 in "uar"?
    ## Note (JH 11/10/18): We should require this.
    uar <- function(p = 2, parcor, sign = "-", atanh.tr = TRUE, fixed = NULL, ...){
        ## Note (JH 11/10/18): adding checks
        ## check p is at least 2
        if(p < 2)
            stop("'p' must be greater than 2 for unit AR polynomials")
        ## the last value should be +/-1
        if(abs(parcor[p]) != 1)
            stop("the last coefficient supplied should be +1 or -1")

        if(is.null(fixed)){
            fixed <- rep(FALSE, p)
            if(p > 0)
                fixed[p] <- TRUE
        }
                            # NOTE: do not name p and parcor in the call for flexibility
        coef <- .set_coef(p, parcor, sign = sign, dispname = "uar", atanh.tr = atanh.tr,
                          fixed = fixed, ...)
                                                
        list(name = "uar", p = p, coef = coef)   # TODO: rename 'coef' to 'parcor' here?
    }

    ## seasonal
    sar <- function(s, p = 0, ar, sign = "-", atanh.tr = TRUE, ...){
                  # NOTE: do not name p and ar in the call for flexibility
        coef <- .set_coef(p, ar, sign = sign, nseasons = s, dispname = "sar",
                          atanh.tr = atanh.tr, ...)
        list(name = "sar", s = s, p = p, coef = coef)        # TODO: remove 's' and 'p'?
    }

    ## Note (JH, 11/10/18): 's' and 'p' aren't used, so removing them won't be a problem.
    sma <- function(s, q = 0, ma, sign = "+", atanh.tr = TRUE, ...){
        coef <- .set_coef(q, ma, sign = sign, nseasons = s, dispname = "sma",
                          atanh.tr = atanh.tr, ...)
        list(name = "sma", s = s, q = q, coef = coef)
    }

    ss <- function(s, ...){
        ## TODO: for now 'm' of length 1
        dots <- c(...)
        if(length(dots) != 1)
            stop("currently only one argument is supported")
        p <- dots[1] - 1
        coef <- .set_coef(p, rep(1, p), sign = "+", nseasons = s, fixed = TRUE,
                          operator = TRUE, dispname = "ss")
        list(name = "ss", s = s, p = p, coef = coef)
    }
    
    ## seasonal unit roots
    si <- function(s, d, ...){
        coef <- .set_coef(1, 1, sign = "-", nseasons = s, fixed = TRUE,
                          operator = TRUE, dispname = "si", ...)
        list(name = "si", s = s, d = d, coef = coef)
    }

    su <- function(s, h){
        # TODO: check that elements of h are among 1,2, ..., s/2 (excluding s/s
        #       if s is even) coef is a list here!
        ## (JH, 11/10/18): Checking h
        int_cond <- any(h %% 1 > 0)
        if(int_cond)
            stop("'h' must contain only (positive) integer values")
        pos_cond <- any(h < 1)
        if(int_cond)
            stop("'h' must contain only positive (integer) values")
        max_cond <- !all(h < s/2)
        if(max_cond)
            stop("'h' must be a positive integer less than 's/2'")

        coef <- lapply(h, function(x)
            .set_coef(2, c(- 2 * cospi(2*x/s), 1),
                      sign = "+", fixed = TRUE, operator = TRUE,
                      dispname = "su", su.nseasons = s, su.harmonic = x) )
        list(name = "su", s = s, u = h, coef = coef) 
    }

    .specials <- ls() # assumes that user facing objects do not start with dot!

    .sarima_descr <- list()
    class(.sarima_descr) <- "sarimadescr"

    environment()
}# end of .sarima_env()

SARIMA <- function(formula){
    e <- .sarima_env()
    parent.env(e) <- environment(formula)
    
                        # keep.order keeps the order for the data frame.
                        # The order in attribute "specials" is not affected by this.
    te <- terms(formula, specials = e$.specials, keep.order = TRUE)

                        # tmp <- (attr(te, "variables")[[3]]); eval(tmp, envir = e)
    termvars <- attr(te, "variables")
    sp <- attr(te, "specials")
    res <- list()

    indices <- lapply(names(sp), 
                      function(nam){ 
                          ind <- sp[[nam]]
                          if(is.null(ind))
                              return(integer(0))
                          names(ind) <- rep(nam, length(ind))
                          ind
                      }
                      )

    indices <- sort(unlist(indices))
    nams <- names(indices)
    
    indices <- indices + 1
    res <- vector("list", length = length(indices))
    names(res) <- nams
    for(i in seq_along(indices)){
        ## evaluate the terms (like ar(), ma(), etc.)
        res[[i]] <- eval(termvars[[indices[i]]], envir = e)
    }

    res
}

trendMaker <- function(formula, data = NULL, time = NULL){
    if(is.null(data))
        data.env <- new.env(hash = FALSE, parent = environment(formula))
    else{ 
        stopifnot(is.environment(data))
        data.env <- data
    }
   
    environment(formula) <- data.env
    data.env$formula <- formula

                 # data.env$.t.orig can be used to detect if 't' is not the original one
    data.env$.t.orig <- data.env$t <- time

    .t.orig <- "dummy for R CMD check; .t.orig is needed in dataenv, as set above"

    data.env$.p <- function(degree){
        if(identical(t, .t.orig))
            res <- poly(t, degree = degree, simple = TRUE) 
        else{
            wrk <- poly(.t.orig, degree = degree, simple = FALSE)
            res <- predict(wrk, newdata = t) 
        }
        colnames(res) <- paste0("ortht", seq_len(degree))
        res
    }
    environment(data.env$.p) <- data.env

        # cs <- function(s, freq){
        #     ## s - scalar; k in 1,2,...,[s/2]
        #     wrk <- lapply(freq, function(k) cbind(cospi(2*time*k/s), sinpi(2*time*k/s)))
        #     res <- do.call("cbind", wrk)
        #     
        #         # colnames(res) <-
        #         #     unlist(lapply(freq, function(x) paste0(c("cos", "sin"), x, ".", s)))
        #     res
        # }
    data.env$.cs <- function(s, k){
        ## s - scalar; k in 1,2,...,[s/2]
        ## TODO: for now assume scalar k too:
        res <- if(length(k) == 1){  # cbind(cospi(2*t*k/s), sinpi(2*t*k/s))
                   structure(c(cospi(2*t*k/s), sinpi(2*t*k/s)), 
                               dim = c(length(t), 2),
                               .Dimnames =  list(NULL, paste0(c("c", "s"), s, ".", k) ))
              }else if(length(k) > 1){
                  wrk <- lapply(k, function(x) cbind(cospi(2*t*x/s), sinpi(2*t*x/s)) )
                  wrk <- do.call("cbind", wrk)
                  colnames(wrk) <- paste0(c("c", "s"), s, ".", rep(k, each = 2) )
                  wrk
              }else
                  stop("'k' must have positive length")

            # TODO: put names, keeping them short (Hm... 7*48 ...)
            # TODO: colnames(res) <- paste0(c("cos", "sin"), k, ".", s)
        res
    }
    environment(data.env$.cs) <- data.env

    ## TODO: test thoroughly
    data.env$.B <- function(x, lags){
        flags <- t > 0
        t.pos <- t[flags]
        t.neg <- t[!flags]

        if(any(flags)){
            ## 15/01/2018 - krapka to accept matrix x
            if(is.matrix(x)){ ## TODO: maybe check also if ncol(1) ? 
                res <- matrix(NA_real_, nrow = length(t), ncol = ncol(x) * length(lags))
                nc.x <- ncol(x)
                curcols <- seq_len(nc.x)
                for(i in seq_along(lags)){
                    lag <- lags[i]
                    res[t - lag > 0, curcols] <- x[pmax(t - lag, 0), ]
                    curcols <- curcols + nc.x 
                }
                ## TODO: don't know why there are no names in this case.
                xnam <- deparse(substitute(x))
                ## TODO: check for "[" or "$" in xnam and adjust the names accordingly
                colnames(res) <- paste0("L(", xnam, rep(seq_len(ncol(x)), length(lags)),
                                        ", ", rep(lags, each = ncol(x)), ")") 
            }else{
                res <- matrix(NA_real_, nrow = length(t), ncol = length(lags))
                for(i in seq_along(lags)){
                    lag <- lags[i]
                    res[t - lag > 0, i] <- x[pmax(t - lag, 0)]
                }
            }
        }
        res
    }
    environment(data.env$.B) <- data.env

    res <- function(index){
        if(missing(index)){# TODO: napravi tova s na.action to avoid losing NA's
            
                # 15/01/2018 was: model.matrix(formula)
                #     changing to use model.frame since we need to keep rows with NA's
                #     and it seems that model.matrix() doesn't have this argument.
                #     TODO: !!! However, model.frame doesn't contain the intercept,
                #               if present! So, it gets lost here.
                #
                # OK, sorted, this is a reply by Brian Ripley to a similar question
                # (TODO: take note also of the user defined function for na.action.):
                #
                #     Yes, but it uses options("na.action"), so you could set that.
                #     Better, call model.frame then model.matrix, something like
                # 
                #     stratmat <- model.matrix(myformula, 
                #                   model.frame(myformula,mydata, na.action = function(x) x))
                # 
                #     ?model.matrix does tell you the second argument should be the result
                #     of model.frame, which is a pretty strong hint.
                 
            res <- model.matrix(formula, model.frame(formula, na.action = NULL))
        }else{
            t.old <- t
            t <<- index
            ## TODO: krapka!!!! - drop the lhs since the length of t is now diferent.

            ## TODO: this will do if there is only time trend but exogeneous vars need
            ##       further work also, will this work if there is only intercept?  Maybe
            ##       should put a dummy variable on the left?
            formula <- formula(Formula::Formula(formula), lhs = FALSE, rhs = 1)
            txt <- paste0(".dummy", paste0(as.character(formula), collapse = ""))
            formula <- as.formula(txt, env = environment(formula))
                # 15/01/2018 - changing to use model.frame since we need to keep rows with
                #     NA's and it seems that model.matrix() doesn't have this argument.
                #     (see also Ripley's remark in  the 'if' clause above)
                # res <- model.matrix(formula, data = data.frame(.dummy = t))
                # res <- as.matrix(res[ , -1])
            res <- model.frame(formula, data = data.frame(.dummy = t), na.action = NULL)
            res <- model.matrix(formula, res)
            t <<- t.old
            res
        }
    }
    environment(res) <- data.env
    res
}

.cat_u_delta <- function(x){
    f <- function(poly){
        s <- environment(poly)$nseasons
        if(!is.null(s)) 
            poly <- poly(poly_x^s)
        if(all(coef(poly) == 1) && (deg <- length(coef(poly)) - 1) > 3)
            res <- paste0("1 + B + ... + B^", deg)
        else
            res <- as.character(poly, "B")
        d <- environment(poly)$d
        if(is.null(d) || d == 1)
            paste0("(", res, ")")
        else
            paste0("(", res, ")^", d)
    }
    paste0(sapply(x, f), collapse = "")
}

.Sarima_fixed <- function(x) x$internal$fixed
.capture_fo <- function(x) if(is.numeric(x)) x  else capture.output(print(x, FALSE))
.Fo.xreg <- function(x) x$internal$Fo.xreg
.Fo.sarima <- function(x) x$internal$Fo.sarima
.Fo.regx <- function(x) if(is.null(res <- x$internal$Fo.regx)) 0 else res

print.Sarima <-
    function (x, digits = max(3L, getOption("digits") - 3L), se = TRUE, ...){
    cat("*Sarima model*")
    cat("\nCall:", deparse(x$call, width.cutoff = 75L), "", sep = "\n")
        # cat("\nTrend: \n\t")
        # print(Fo.trend, FALSE)
        # cat("\nSARIMA specification: \n\t")
        # print(Fo.sarima, FALSE)
    delta <- x$internal$delta
    if(length(delta) > 0){
        cat("Unit root terms:\n", "    ")
        cat(.cat_u_delta(x$internal$delta_poly), sep = "")
        cat("\n\n")
    }

    if (length(x$coef)) {
        cat("Coefficients:\n")
        coef <- round(x$coef, digits = digits)
        ## use NROW as if all coefs are fixed there are no var.coef's
        if (se && NROW(x$var.coef)) {
            ses <- rep.int(0, length(coef))
            ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
            coef <- matrix(coef, 1L, dimnames = list(NULL, names(coef)))
            coef <- rbind(coef, s.e. = ses)
        }
        print.default(coef, print.gap = 2)
    }
    cm <- x$call$method
    if(is.null(cm) || cm != "css")
        cat("\nsigma^2 estimated as ", format(x$sigma2, digits = digits),
            ":  log likelihood = ", format(round(x$loglik, 2L)),
            ",  aic = ", format(round(x$aic, 2L)), "\n", sep = "")
    else
        cat("\nsigma^2 estimated as ",
            format(x$sigma2, digits = digits),
            ":  part log likelihood = ", format(round(x$loglik,2)),
            "\n", sep = "")
    invisible(x)
}

summary.Sarima <- function(object, ...){
    ## follow the 'lm' tradition
    cat("\nCall:\n", paste(deparse(object$call), sep = "\n", collapse = "\n"), 
        "\n\n", sep = "")

    ## TODO: these need more stable implementation, probably should be in the object.
               # object$call[[2]][[3]][[1]] is '|'
        # Fo.xreg   <- object$formula[[3]][[2]] # or object$call[[2]][[3]][[2]]
        # Fo.sarima <- object$formula[[3]][[3]]
        # Fo.regx   <- if(length(object$formula[[3]]) > 3) object$formula[[3]][[4]] else 0

    Fo.xreg   <- .Fo.xreg  (object)
    Fo.sarima <- .Fo.sarima(object)
    Fo.regx   <- .Fo.regx  (object)

    xfo      <- .capture_fo(Fo.xreg)
    fosarima <- .capture_fo(Fo.sarima)
    fox      <- .capture_fo(Fo.regx)
      
    cat("Model:  Y_t - xreg_t is SarimaX\n")
    cat("  xreg:   ", xfo, "\n")
    cat("  sarima: ", fosarima, "\n")
    cat("  regx:   ", fox, "\n")
    cat("\n")

    delta <- object$internal$delta
    if(length(delta) > 0){
        cat("Unit root terms:\n", "    ")
        cat(.cat_u_delta(object$internal$delta_poly), sep = "")
        cat("\n\n")
    }

    fixed <- .Sarima_fixed(object)
    est <- object$coef[!fixed]
    se <- sqrt(diag(object$var.coef))[!fixed]
    zval <- est / se
    pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
    coefs <- cbind(est, se, zval, pval)
    dimnames(coefs) <- 
        list(names(est), 
             c("Estimate", "Std. Error", "Z value", "Pr(>|z|)") )

    cat("Coefficients:\n")
    printCoefmat(coefs
                 # , digits = digits, signif.stars = signif.stars, na.print = "NA", ...
                )

    cat("\nestimated sigma^2 = ", format(object$sigma2, digits = 3), 
        ",  log-likelihood = ", format(round(object$loglik, 2)),       
        ",  aic = ", format(round(object$aic, 2L)),
        sep = "")

#browser()
    cat("\n")

    invisible(object)
}

.tsdiag_choices <- c(
    ## "classic (std. residuals, acf, portmanteau p-values)",
    "residuals",
    "acf of residuals",
    "p values for Ljung-Box statistic",
    "p values for Li-McLeod statistic",
    "p values for Box-Pierce statistic",
    "pacf of residuals"
)

tsdiag.Sarima <- function(object, gof.lag = NULL, ask = FALSE, ..., plot = 1:3, layout = NULL)
{
    if(is.null(gof.lag))
        gof.lag <- 20  # NOTE: arbitrary value
    else if(!is.numeric(gof.lag))
        stop("'gof.lag' must be numeric and contain positive integers")

    lag.max <- max(gof.lag)


    ## plot standardized residuals, acf of residuals, Ljung-Box p-values
    err <- object$residuals
    stdres <- err / sqrt(object$sigma2)

    choices <- .tsdiag_choices
    chnum <- 1:length(choices)
    
    if(!isTRUE(plot)){                  # plot is typically numeric index here;
        choices <- choices[plot]        # FALSE or NULL give zero length result, so no plots
        chnum <- chnum[plot]

        if(anyNA(choices)){
            warning("'plot' should be TRUE/FALSE or vector of positive integers <= ",
                    length(.tsdiag_choices), ",\n", "ignoring non-existent values")
            chnum <- chnum[!is.na(choices)]
            choices <- choices[!is.na(choices)]
        }
    }

    if(length(choices) > 0){
	old.par <- par(no.readonly = TRUE)
	on.exit(par(old.par))     # restore graphics parameters before exiting.

        # n_per_page <- 3

        n_per_page <- if(is.null(layout))
                          layout(matrix(1:min(3, length(choices)), ncol = 1))
                      else
                          ## TODO: this needs further thought!
                          do.call("layout", layout) # for the time being

        ask_user <- interactive() && (ask || length(choices) > n_per_page)
            
        choice_title <- "Select a plot number or 0 to exit"
        ch_index <- if(length(choices) == 1)
                        1
                    else if(ask)
                        menu(choices, title = choice_title)
                    else if(!identical(plot, FALSE))
                        1
                    else
                        integer(0)
        choice <- chnum[ch_index]

        ## precompute common stuff for portmanteau tests
        nlag <- gof.lag
        pval <- numeric(nlag)
        fitdf <- if(inherits(object, "Sarima"))
                     length(object$internal$nonfixed)
                 else if(inherits(object, "Arima"))
                     sum(object$arma[1:4]) # object$arma is: p, q, p_s. q_s, s, d, d_s
                 else
                     0
                    # for(i in 1L:nlag)
                    #     pval[i] <- Box.test(err, i, type="Ljung-Box", 
                    #                         fitdf = ifelse(i > fitdf, fitdf, i - 1))$p.value
        sacf <- autocorrelations(err, maxlag = nlag) # deal with NA's?
        
        res <- list(residuals = err,
                    "LjungBox"   = NULL,
                    "LiMcLeod"   = NULL,
                    "BoxPierce" = NULL)
        while(length(choice) != 0){
            switch(choice,
            { # 1:  "residuals",
                plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
                abline(h = 0)
                #acf(err, main = "ACF of residuals from model", lag.max = lag.max)
                #pacf(err, main = "PACF of residuals from model", lag.max = lag.max)
            },
            { # 2:  "ACF of residuals"
                ## acf
                acf(err, plot = TRUE, main = "ACF of Residuals", lag.max = lag.max, na.action = na.pass)
                    # acf(cdf, main = "", lag.max = lag.max)
                    # title("ACF of" ~U[t])
                    # pacf(cdf, main = "", lag.max = lag.max)
                    # title("PACF of" ~U[t])
            },
            { # 3: "Ljung-Box p-values"
		acftest <- acfIidTest(sacf, npar = fitdf, nlags = 1:nlag, method = "LjungBox",
                                      interval = NULL)
                res[["LjungBox"]] <- acftest
            },
            { # 4: "Li-McLeod p-values"
		acftest <- acfIidTest(sacf, npar = fitdf, nlags = 1:nlag, method = "LiMcLeod",
                                      interval = NULL)
                res[["LiMcLeod"]] <- acftest
            },
            { # 5: "Box-Pierce p-values"
		acftest <- acfIidTest(sacf, npar = fitdf, nlags = 1:nlag, method = "BoxPierce",
                                      interval = NULL)
                res[["BoxPierce"]] <- acftest
            },
            { # 6:  "PACF of residuals"
                ## acf
                pacf(err, plot = TRUE, main = "PACF of Residuals", lag.max = lag.max, na.action = na.pass)
            },
                # { # 4: "ACF/Histogram of tau_residuals"
                #     acf(err2, main = "ACF of tau_residuals", lag.max = lag.max)
                #     hist(err2, freq = FALSE, main = "Histogram of tau_residuals", xlab  =  "",
                #          ylim = c(0, 0.5))
                #     lines(seq(-5, 5, .01), dnorm(seq(-5, 5, .01)), col = "red")
                # }
            )

            if(choice %in% 3:5){ # common code for portmanteau tests
                pval <- acftest$test[ , "pvalue"]
                plot(1L:nlag, pval, xlab = "lag", ylab = "p value", ylim = c(0,1),
                     main = .tsdiag_choices[choice])
                abline(h = 0.05, lty = 2, col = "blue")
            }
            
            if(length(chnum) == 1)  # length(choices) == 1
                break

            if(ask_user) {
                ch_index <- menu(choices, title = choice_title)
                choice <- chnum[ch_index]
            } else{
                ## just plot the next one
                ##  Note: this doesn't update ch_index
                chnum <- chnum[-1]
                choice <- chnum[1]
            }
        }
    }

    class(res) <- "tsdiagSarima"
    
    invisible(res)
}

## fkf: y ~ d x 1
##      state ~ m x 1
armapqss <- function(ar, ma, sigma) {
    p <- length(ar)
    q <- length(ma)
    r <- max(p, q + 1) # fkf: m

    ear <- c(ar, numeric(max(r - p, 0)))
    ema <- c(ma, numeric(max(r - q - 1, 0)))

    Tt <- cbind(ear, rbind(diag(r - 1), 0))
    Rt <- matrix(c(1, ema), ncol = 1)

    Zt <- matrix(c(1, numeric(r - 1)), nrow = 1)

    ct <- matrix(0)                      # obs. regression
    dt <- matrix(0, nrow = r, ncol = 1)  # state regression

    GGt <- matrix(0)    # d x d

    H <- Rt * sigma
    HHt <- H %*% t(H)

    a0 <- numeric(r)
    P0 <- diag(1e6, nrow = r)
    
    list(a0 = a0, P0 = P0, ct = ct, dt = dt,
         Zt = Zt, Tt = Tt, GGt = GGt, HHt = HHt)
}

## extention of armapqss
## fkf: y ~ d x 1
##      state ~ m x 1
xarmaxss <- function(ar, ma, sigma, xreg, regx) {
    d <- 1  # dim of y, currently univariate

    p <- length(ar)
    q <- length(ma)
    r <- max(p, q + 1) # fkf: m = r

    ear <- c(ar, numeric(max(r - p, 0)))
    ema <- c(ma, numeric(max(r - q - 1, 0)))

    Tt <- cbind(ear, rbind(diag(r - 1), 0))
    Rt <- matrix(c(1, ema), ncol = 1)

    Zt <- matrix(c(1, numeric(r - 1)), nrow = d)

    ## TODO: more flexible options for ct i dt ?

    ## obs. regression
    ct <- if(missing(xreg) || length(xreg) == 0) 
              matrix(0, nrow = d)
          else
              xreg

    ## state regression
    dt <- if(missing(regx) || length(regx) == 0)
              matrix(0, nrow = r, ncol = 1)
          else  if(nrow(regx) == r)
              regx
          else # assume < r
              rbind(regx, matrix(0, nrow = r - nrow(regx), ncol = ncol(regx)))

    GGt <- matrix(0, nrow = d)    # d x d

    H <- Rt * sigma
    HHt <- H %*% t(H)

    a0 <- numeric(r)
    P0 <- diag(1e6, nrow = r)
    
    list(a0 = a0, P0 = P0, ct = ct, dt = dt,
         Zt = Zt, Tt = Tt, GGt = GGt, HHt = HHt)
}

## extention of xarmaxss
## fkf: y ~ d x 1
##      state ~ m x 1
xarimaxss <- function(ar, ma, sigma, xreg, regx, delta = numeric(0)) {
    model <- xarmaxss(ar, ma, sigma, xreg, regx)
    dd <- length(delta)
    if(dd == 0)
        return(model)
    
    model$dt <- c(model$dt, numeric(dd)) # extend dt with zeroes ## TODO: dt is a matrix!
    model$Tt <- rbind(c(delta, 1, numeric(ncol(model$Tt) - 1)),
                      dbind(diag(1, nrow = dd - 1, ncol = dd),
                            model$Tt ))
    model$Ht <- dbind(matrix(0, dd, dd), model$Ht)
    model$Zt <- c(model$Zt, numeric(dd))   # [1, 0, ...., 0] 

    model$a0 <- c(numeric(dd), model$a0)
    model$P0 <- dbind(diag(1e6, nrow = dd), model$P0)
    
        # list(a0 = a0, P0 = P0, ct = ct, dt = dt,
        #      Zt = Zt, Tt = Tt, GGt = GGt, HHt = HHt)
    model
}

## modification of makeARIMA() from ~R/src/base/R-3.3.2/src/library/stats/R/~
makeArimaGnb <- function(phi, theta, Delta, kappa = 1e6,
                      SSinit = "gnb", # c("Rossignol2011", "gnb", "Gardner1980")
                      tol = .Machine$double.eps){
    if(anyNA(phi))   warning(gettextf("NAs in '%s'", "phi"), domain=NA)
    if(anyNA(theta)) warning(gettextf("NAs in '%s'", "theta"), domain=NA)
    p <- length(phi)
    q <- length(theta)
    d <- length(Delta)

    r <- max(p, q + 1L)
    rd <- r + d
    V <- Pn <- P <- T <- matrix(0., rd, rd)
    h <- 0.
    a <- rep(0., rd)
    Z <- c(1., rep.int(0, r-1L), Delta)
    
    ## if(q < r - 1L)
    ##     theta <- c(theta, rep.int(0, r-1L-q))
    ## R <- c(1, theta, rep.int(0, d))
    ## V <- R %o% R
    if(q > 0)
        V[1:(q+1), 1:(q+1)] <- c(1, theta) %o% c(1, theta)
    else  
        V[1, 1] <- 1.0

    if(q < r - 1L)  # extent theta - this is in the returned value!
        theta <- c(theta, rep.int(0, r - 1L - q))
    
    ## The stationary part of T is left-companion of dim (r,r)
    if(p > 0)
        T[1L:p, 1L] <- phi
    if(r > 1L) {
        ind <- 2:r
        T[cbind(ind-1L, ind)] <- 1
    }
    ## set the top left corner of Pn to the covariance matrix of the stationary part
    if(r > 1L)
        Pn[1L:r, 1L:r] <- switch(match.arg(SSinit),
                                 # "Gardner1980" = .Call(stats:::C_getQ0, phi, theta),
                                 # "Rossignol2011" = .Call(stats:::C_getQ0bis, phi, theta, tol),
                                 "gnb" = .Call(`_sarima_arma_Q0gnb`, phi, theta, tol),
                                 stop("invalid 'SSinit'"))
    else Pn[1L, 1L] <- if(p > 0) 1/(1 - phi^2) else 1

    if(d > 0L) {
        ## The nonstationary part of T is top-companion of dim (d,d);
        ##      the r+1 row of T is equal to Z: its first element is 1,
        ##      followed by zeroes, followed by Delta from column r + 1.
        T[r + 1L, ] <- Z
        if(d > 1L) {
            ind <- r + 2:d
            T[cbind(ind, ind-1)] <- 1
        }
        ## set the diagonal of the bottom right corner of Pn to kappa
        ##     (this corresponds to the nonstationary part)
        Pn[cbind(r + 1L:d, r + 1L:d)] <- kappa
    }

    ## note P is left equal to the zero matrix
    list(phi = phi, theta = theta, Delta = Delta, Z = Z, a = a,
         P = P, T = T, V = V, h = h, Pn = Pn)
}

sarimaReport <- function(o1, o2, ...){

    list(
        all.equal(coef(o1), coef(o2))
        )


}

factorizeMA <- function(x, theta = NULL, tol = 1e-12, maxiter = 1000){
    q <- length(x) - 1
    if(is.null(theta))
        theta <- c(sqrt(x[1] + 2 * sum(x[-1])), rep(0, q))

    storage.mode(theta) <- "double"

        # r <- ltsa::tacvfARMA(theta = -theta[-1]/theta[1], maxLag = q, sigma2 = theta[1]^2)
    r <- .Call("_sarima_MAacvf0", theta)
#cat( "theta = ", theta, ",\t r = ", r, "\n")

    crit <- 1e100
    iter <- 0
    while(crit > tol  &&  iter <= maxiter){
        iter <- iter + 1
        T <- .Call("_sarima_DAcvfWrtMA", theta)

        theta <- solve(T, r + x)

        r <- ltsa::tacvfARMA(theta = -theta[-1]/theta[1], maxLag = q, sigma2 = theta[1]^2)

        rnew <- .Call("_sarima_MAacvf0", as.double(theta)) # MAacvf0(theta)

#cat( "theta = ", theta, ",\t r = ", r, "\n")
#browser()
        crit <- sum((x - r)^2)
    }

    ## output as SQUAREM::fpiter(), maybe should use it anyway.
    list(par = theta, value.objfn = crit, fpevals = iter, 
        objfevals = iter, convergence = if(iter > maxiter) 1 else 0)

}

##  arimaSS <- function(y, mod)
##     {
##         ## next call changes mod components a, P, Pn so beware!
##         .Call(stats:::C_ARIMA_Like, y, mod, 0L, TRUE)
##     }