1
$\begingroup$

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?

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

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"))
}
$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Frank Feb 1 '12 at 16:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.