[R-sig-finance] tips and tricks for rolling regressions?

[Ajay Shah]

On Wed, Aug 25, 2004 at 12:31:51PM +0100, Taher Khan wrote:
> Greetings!
> 
> I want to find out if there is a clever way to do a one-step-ahead
> moving window rolling regression using a vectorized operation like with
> apply, rather than a for loop.
> 
> Does anyone have a code snippet they would like to share for how they do
> rolling regressions?!?!

Version 1 : my low-quality hacked-up version -- but using simple
programming constructs that even a child can understand --

snip snip ------------------------------------------------------------
# Using moving windows, we do an OLS regression.
#
# The function returns a matrix with T rows. So there are plenty of
# "NA" values in it. The columns are organised as:
#    beta1 s1 beta2 s2 beta3 s3 sigma R^2
# i.e. we have each regression coefficient followed by it's sigma,
# and then at the end we have the residual sigma and the regression R^2.
#
# Example of usage -- to get moving window mean & sigma - 
#  r = movingWindowRegression(E, length(E$date), width, r.usd ~ 1, 1)
# WARNING: for the model, say 'r.usd~1', not 'E$r.usd~1'.
#
# Use the returned matrix like r[width:T,1] as the moving window mean etc.
movingWindowRegression <- function(data, T, width, model, K) {
  results = matrix(nrow=T, ncol=2*K+2)
  for (i in width:T) {
    details <- summary.lm(lm(model, data[(i-width+1):i,]))
    n=1;
    for (j in 1:K) {
      results[i, n:(n+1)]   = details$coefficients[j, 1:2]
      n = n + 2
    }
    results[i, n] = details$sigma
    results[i, n+1] = details$r.squared
  }
  return(results)
}
snip snip ------------------------------------------------------------






And from Douglas Bates (bates at stat.wisc.edu), a beautiful, brilliant
version, that is faster than mine, but that I only dimly understand --
see http://tolstoy.newcastle.edu.au/R/help/04/04/1269.html -- 

snip snip ------------------------------------------------------------
movingWindow2 <- function(formula, data, width, ...) {
    mCall = match.call()
    mCall$width = NULL
    mCall[[1]] = as.name("lm")
    mCall$x = mCall$y = TRUE
    bigfit = eval(mCall, parent.frame())
    ncoef = length(coef(bigfit))
    width = as.integer(width)[1]
    y = bigfit$y
    x = bigfit$x
    nr = nrow(x)
    stopifnot(width >= ncoef, width <= nr)
    terms = bigfit$terms
    inds = embed(seq(nr), width)[, rev(seq(width))]
    sumrys <- lapply(seq(nrow(inds)),
                     function(st) {
                         ind = inds[st,]
                         fit = lm.fit(x[ind, , drop = FALSE], y[ind])
                         fit$terms = terms
                         class(fit) = "lm"
                         summary(fit)
                     })
    list(coefficients = sapply(sumrys, function(sm) coef(sm)[,"Estimate"]),
         Std.Error = sapply(sumrys, function(sm) coef(sm)[,"Std. Error"]),
         sigma = sapply(sumrys, "[[", "sigma"),
         r.squared = sapply(sumrys, "[[", "r.squared"))
}
snip snip ------------------------------------------------------------

-- 
Ajay Shah                                                   Consultant
ajayshah at mayin.org                      Department of Economic Affairs
http://www.mayin.org/ajayshah           Ministry of Finance, New Delhi