Hi This is my first post to this group, so apologies in advance if I get it wrong. I would like to know how the prediction for arima models works in R. I have a time series to which I fit an arima model, of varying AR and MA orders. I then use the predict function to project it forward. I have also written my own function to perform the prediction, but it gives different answers to Arima.predict when the MA order is non-zero. I use the residuals from the arima function in my custom prediction function. I think this may be my problem. In the arima model: x{t} = a(1)x{t-1} + a(2)x{t-2} + ... + a(p)x{t-p} + e{t} + b(1)e{t-1} + b(2)e{t-2} + ... + b(q)e{t-q} I am treating the residuals (i.e. arima(....)$res) as the e{t} terms. This gives different answers both in the region of the simulation and in the region of the prediction, so I'm guessing they're not what I think they are. Indeed, after q intervals in the prediction, the prediction proceeds as I would expect, presumably because all the residuals that have an effect are zero by this stage. Any help greatly appreciated - my code is below. Thanks Chris -------------------------------------------------------------- The code to produce the two predictions is as follows: AR <- 5 MA <- 3 sim <- arima.sim(list(order=c(AR,0,MA), ar=c(.1, .1, .1, .1, .1), ma=c(.1, .1, .1)), n=100) + 50 fit <- arima(sim, order = c(AR, 0, MA)) coefs <- fit$coef series <- sim innov <- fit$res pred <- 100 fit.predict <- predict(fit, n.ahead = pred) fit.r <- c(sim, fit.predict$pred) fit.custom <- ProjectCentralArima(AR = AR, MA = MA, d = 0, coefs coefs, series = sim, innov = innov, pred = pred)$ser ProjectCentralArima function: ProjectCentralArima <- function(AR, MA, d, coefs, series, innov, pred) { if(d==0){ series.diff <- series } else { series.diff <- diff(series, lag=1, differences=d) } intercept <- coefs[length(coefs)] for(i in 1:pred){ temp <- intercept l.s <- length(series.diff) if(AR > 0){ for(j in 1:AR){ temp <- temp + coefs[j] * (series.diff[l.s - j + 1] - intercept) } } if(MA > 0){ for(j in (1:MA)){ temp <- temp + coefs[j + AR] * innov[l.s - j + 1] } } innov <- c(innov, 0) series.diff <- c(series.diff, temp) } if(d==0){ series.undiff <- series.diff } else { series.undiff <- diffinv(series.diff, lag=1, differences=d, xi series[1:d]) } return(list(series = series.undiff, innov = innov)) }