# Random walk out of sample forecasting

Probably my question is a bit stupid, but I'm having some problems in writing down in R the out-of-sample forecasting with a Random Walk. I have a multivariate time series (y) and I want to estimate the out of sample forecasting (y(t+k)-y^(t+k/t) result with a RW, for k =1,6,12. I write down this code, it works, but it seems to give me too much low errors.

residuals1 <- rep(0,58)
residuals6 <- rep(0,58)
residuals12 <- rep(0,58)

y1 <- t(y[,1])
for (i in 1:58) {
residuals1[i] <- y1[134+i+1]-y1[134+i]
residuals6[i] <- y1[134+i+6]-y1[134+i]
residuals12[i] <- y1[134+i+12]-y1[134+i]
}


I think that I'm doing some mistakes, do you have any suggestions?

• You may want to rewrite the expression for the out-of-sample forecast, possibly using LaTeX. – chl Feb 1 '12 at 12:32

The rwf function in the forecast package may be useful. You can either use it directly, or take a look at the source code:

> library(forecast)
> rwf
function (x, h = 10, drift = FALSE, level = c(80, 95), fan = FALSE,
lambda = NULL)
{
xname <- deparse(substitute(x))
n <- length(x)
freq = frequency(x)
nn <- 1:h
if (!is.ts(x))
x <- ts(x)
if (!is.null(lambda)) {
origx <- x
x <- BoxCox(x, lambda)
}
if (drift) {
fit <- summary(lm(diff(x) ~ 1))
b <- fit$coefficients[1, 1] b.se <- fit$coefficients[1, 2]
s <- fit$sigma res <- ts(c(NA, residuals(fit))) method <- "Random walk with drift" } else { b <- b.se <- 0 s <- sd(diff(x), na.rm = TRUE) res <- ts(c(NA, diff(x))) method <- "Random walk" } tsp(res) <- tsp(x) f <- x[n] + nn * b se <- sqrt((nn * s^2) + (nn * b.se)^2) if (fan) level <- seq(51, 99, by = 3) else { if (min(level) > 0 & max(level) < 1) level <- 100 * level else if (min(level) < 0 | max(level) > 99.99) stop("Confidence limit out of range") } nconf <- length(level) z <- qnorm(0.5 + level/200) lower <- upper <- matrix(NA, nrow = h, ncol = nconf) for (i in 1:nconf) { lower[, i] <- f - z[i] * se upper[, i] <- f + z[i] * se } lower <- ts(lower, start = tsp(x)[2] + 1/freq, f = freq) upper <- ts(upper, start = tsp(x)[2] + 1/freq, f = freq) colnames(lower) <- colnames(upper) <- paste(level, "%", sep = "") fcast <- ts(f, start = tsp(x)[2] + 1/freq, f = freq) fits <- x - res if (!is.null(lambda)) { x <- origx fcast <- InvBoxCox(fcast, lambda) fits <- InvBoxCox(fits, lambda) upper <- InvBoxCox(upper, lambda) lower <- InvBoxCox(lower, lambda) } junk <- list(method = method, level = level, x = x, xname = xname, mean = fcast, lower = lower, upper = upper, model = list(drift = b, drift.se = b.se, sd = s), fitted = fits, residuals = res, lambda = lambda) junk$model\$call <- match.call()
return(structure(junk, class = "forecast"))
}

• Thanks. your answer probably exceeds my needs cause I'm just trying to estimate the out of sample forecasting errors, and using a RW model, the forecast k step ahead y(t+k) it's equal to y(t).And this is what I wrote in my code but I obtained too much small errors, mach smaller than the paper I'm trying to replicate. For this reason I think that I did a mistake. But I can't find it! – Frank Feb 1 '12 at 16:29