Random walk out of sample forecasting

Question

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?

Answer 1

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"))
}