R/predict.Sarima.R

In sarima: Simulation and Prediction with Seasonal ARIMA Models

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

## see predict.arima also
predict.Sarima <- function(object, n.ahead = 1, newxreg = NULL, 
                           se.fit = TRUE, newdata = NULL, ...){
    # myNCOL <- function(x) if (is.null(x)) 
    #     0
    # else NCOL(x)
    # rsd <- object$residuals
    # xr <- object$call$xreg
    # xreg <- if (!is.null(xr)) 
    #     eval.parent(xr)
    # else NULL
    # ncxreg <- myNCOL(xreg)
    # if (myNCOL(newxreg) != ncxreg) 
    #     stop("'xreg' and 'newxreg' have different numbers of columns")
    # class(xreg) <- NULL
    # xtsp <- tsp(rsd)
    # n <- length(rsd)
    # arma <- object$arma
    # coefs <- object$coef
    # narma <- sum(arma[1L:4L])
    # if (length(coefs) > narma) {
    #     if (names(coefs)[narma + 1L] == "intercept") {
    #         xreg <- cbind(intercept = rep(1, n), xreg)
    #         newxreg <- cbind(intercept = rep(1, n.ahead), newxreg)
    #         ncxreg <- ncxreg + 1L
    #     }
    #     xm <- if (narma == 0) 
    #         drop(as.matrix(newxreg) %*% coefs)
    #     else drop(as.matrix(newxreg) %*% coefs[-(1L:narma)])
    # }
    # else 
    ## TODO: no xreg yet
    n <- length(object$series)
    nco <- object$internal$n.coef
    # co.ind <- cumsum(c(0, nco))
    # 2018-02-17 repalce numeric indices in nco[] with character
    arma.ind <- seq_len(nco["arma"])
    xreg.ind <- nco["arma"] + seq_len(nco["xreg"])
    regx.ind <- nco["arma"] + nco["xreg"] + seq_len(nco["regx"])
    if(length(xreg.ind) == 0){
        xm <- 0
    }else{
        trmake <- object$internal$trendMaker
        xregmat <- trmake(n + 1:n.ahead)
        xreg.coef <- object$coef[xreg.ind]
        xm <- as.vector(xregmat %*% xreg.coef)
    }
    ## TODO: no regx yet
    # if (arma[2L] > 0L) {
    #     ma <- coefs[arma[1L] + 1L:arma[2L]]
    #     if (any(Mod(polyroot(c(1, ma))) < 1)) 
    #         warning("MA part of model is not invertible")
    # }
    # if (arma[4L] > 0L) {
    #     ma <- coefs[sum(arma[1L:3L]) + 1L:arma[4L]]
    #     if (any(Mod(polyroot(c(1, ma))) < 1)) 
    #         warning("seasonal MA part of model is not invertible")
    # }
    
    if(is.null(newdata)){
        #--- Can be done original way if no newdata
        z <- sarima_KalmanForecast(n.ahead, object$model)
        pred <- z[[1L]] + xm
        if (se.fit) {
            se <- sqrt(z[[2L]] * object$sigma2)
            pred <- list(pred = pred, se = se)
        }
    }else{
        #
        #z <- NULL
        #mod <- object$model
        #for(i in seq_along(newdata)){
        #    zz <- KalmanForecast(n.ahead, mod)
        #    #--- This only returns the desired n.step ahead - do we want more?
        #    z$pred <- c(z$pred, zz$pred[n.ahead])
        #    z$var <- c(z$var, zz$var[n.ahead]) 
        #    z0 <- KalmanLike(newdata[i], mod, nit = 0L, update = TRUE)
        #    mod <- attr(z0, "mod")
        #}
        #
        z <- sarima_KalmanForeUp(newdata, n.ahead, object$model)
        pred <- lapply(z[[1L]], function(x) x + xm)
        if (se.fit) {
            se <- lapply(z[[1L]], function(x) sqrt(x * object$sigma2))
            pred <- list(pred = pred, se = se)
        }
    }
    ## NOTE (JH, 15/10/18): This assumes no lagged x-values are used! 
    #--- Returns the updated model too
    #--- If original method used, this will be NULL
    pred
}