[R] What's the predict procedure of ARIMA in R?

[Saji Ren]

Hello,guys:
Recently, I am working on a seasonal ARIMA model. And I met some problem in the forecasting.
Now I just want to know that How does R perform the predict procedure(the predict formula, the initial setting of errors,etc.)?
I run the following commands and get the original code of the "predict" command, but I can't read it.
Can anybody explain it to me?

Thanks!
saji from Shanghai

the code:
> getS3method("predict","Arima")
function (object, n.ahead = 1, newxreg = NULL, se.fit = TRUE, 
    ...) 
{
    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[1:4])
    if (length(coefs) > narma) {
        if (names(coefs)[narma + 1] == "intercept") {
            xreg <- cbind(intercept = rep(1, n), xreg)
            newxreg <- cbind(intercept = rep(1, n.ahead), newxreg)
            ncxreg <- ncxreg + 1
        }
        xm <- if (narma == 0) 
            drop(as.matrix(newxreg) %*% coefs)
        else drop(as.matrix(newxreg) %*% coefs[-(1:narma)])
    }
    else xm <- 0
    if (arma[2] > 0) {
        ma <- coefs[arma[1] + 1:arma[2]]
        if (any(Mod(polyroot(c(1, ma))) < 1)) 
            warning("MA part of model is not invertible")
    }
    if (arma[4] > 0) {
        ma <- coefs[sum(arma[1:3]) + 1:arma[4]]
        if (any(Mod(polyroot(c(1, ma))) < 1)) 
            warning("seasonal MA part of model is not invertible")
    }
    z <- KalmanForecast(n.ahead, object$model)
    pred <- ts(z[[1]] + xm, start = xtsp[2] + deltat(rsd), frequency = xtsp[3])
    if (se.fit) {
        se <- ts(sqrt(z[[2]] * object$sigma2), start = xtsp[2] + 
            deltat(rsd), frequency = xtsp[3])
        return(list(pred = pred, se = se))
    }
    else return(pred)
}
<environment: namespace:stats>