Эффективный способ создания рекурсивных оценок вне выборки для расчета RMSE в R

Question

Я пытаюсь оценить производительность прогнозирования вне выборки для разных моделей OLS. Самый простой временной ряд регрессии выглядит следующим образом: Y_t = b0 + b1 * Y_t-30 + e_t

Фитинг период для модели, скажем, 50, то я позволил модели запустить с помощью пакета dynlm

dynlm(as.zoo(Y) ~ L(as.zoo(Y), 30), start = "1996-01-01", end = timelist[i]) 

В моем текущем коде я просто разрешаю индексировать i до конца, а затем сохраняю RMSE соответствующей модели. Но это RMSE не является предварительным прогнозом на один шаг вперед, и поскольку мой текущий код уже довольно медленный, и он даже не делает то, что я хочу, я хотел бы спросить вас, есть ли у вас предложение, которое пакет, который я должен использовать для достижения своей цели.

Подводя итог, я хочу сделать следующее:

1) запустить рекурсивные регрессии после определенного фитинга периода (расширение окна, не прокатит окно)

2) создать один шаг вперед вне образца прогнозы

3) вычислить корень среднеквадратичной погрешности этих прогнозов против фактических наблюдений для оценки эффективности модели

Я пытался делать это до сих пор с огромным для цикла и пакет dynlm , но res Ульты не очень удовлетворяют. Любой вход очень ценится, так как я давно искал решения. Я буду обновлять свой примерный код, как только я достиг некоторого прогресса.

# minimal working example 
require(xts) 
require(zoo) 
require(dynlm) 
timelist <- seq.Date(from = as.Date("1996-01-01"), to = as.Date("1998-01-01"), by = "days") 
set.seed(123) 
Y <- xts(rnorm(n = length(timelist)), order.by = timelist) 
X <- xts(rnorm(n = length(timelist), mean = 10), order.by = timelist) 
# rmse container 
rmse.container.full <- data.frame(matrix(NA, ncol = 3, nrow = length(index(timelist)))) 
colnames(rmse.container.full) <- c("Date", "i", "rmse.m1") 
rmse.container.full$Date <- timelist 
# fitting period 
for(i in 50:length(timelist)) { 
    # m1 
    model1 <- dynlm(as.zoo(Y) ~ L(as.zoo(X), 30), start = "1996-01-01", end = timelist[i]) 
    rmse.container.full[i, 2] <- i 
    rmse.container.full[i, 3] <- summary(model1)$sigma # RSME mod1 etc 
    print(i) 
} 

Answer 1

Ну, как тот, кто задал этот вопрос, я хотел бы внести свой вклад, как я решить мою проблему:

Как мне нужно только один шаг вперед прогнозы я могу выбросить все остальное, и это сделало быстрее запустить код. (от 12 минут до ~ 10 секунд на модель).

Я сам создал полный кадр (включая лага) и использовал lm вместо dynlm. Следующий код дал мне мои желаемые результаты (я проверил первые несколько наблюдений вручную и, похоже, работает). Код адаптирован здесь: Recursive regression in R

 mod1.predictions <- lapply(seq(1400, nrow(df.full)-1), function(x) { 
       mod1 <- lm(Y ~ X, data = df.full[1:x, ]) 
       pred1 <- predict(mod1, newdata = df.full[x+1, ]) 
       return(pred1) 
       }) 

Для вычисления RMSE я использовал эту функцию

# rmse function 
rmse <- function(sim, obs) { 
    res <- sqrt(mean((sim - obs)^2, na.rm = TRUE)) 
    res 
} 

Спасибо за намеки на CrossValidated, это помогло много.

Answer 2

Вы можете уменьшить время вычислений достаточно много, используя Fortran функции из

Миллер, А. J. ​​(1992). Алгоритм AS 274: Процедуры с наименьшими квадратами до дополняют функции Джентльмена. Journal of the Royal Statistics Общество. Серия C (прикладная статистика), 41 (2), 458-478.

Вы можете сделать это, используя этот код

# simulate data 
set.seed(101) 
n <- 1000 
X <- matrix(rnorm(10 * n), n, 10) 
y <- drop(10 + X %*% runif(10)) + rnorm(n) 
dat <- data.frame(y = y, X) 

# assign wrapper for biglm 
biglm_wrapper <- function(formula, data, min_window_size){ 
    mf <- model.frame(formula, data) 
    X <- model.matrix(terms(mf), mf) 
    y - model.response(mf) 

    n <- nrow(X) 
    p <- ncol(X) 
    storage.mode(X) <- "double" 
    storage.mode(y) <- "double" 
    w <- 1 
    qr <- list(
    d = numeric(p), rbar = numeric(choose(p, 2)), 
    thetab = numeric(p), sserr = 0, checked = FALSE, tol = numeric(p)) 
    nrbar = length(qr$rbar) 
    beta. <- numeric(p) 

    out <- numeric(n - min_window_size - 2) 
    for(i in 1:(n - 1)){ 
    row <- X[i, ] # will be over written 
    qr[c("d", "rbar", "thetab", "sserr")] <- .Fortran(
     "INCLUD", np = p, nrbar = nrbar, weight = w, xrow = row, yelem = y[i], 
     d = qr$d, rbar = qr$rbar, thetab = qr$thetab, sserr = qr$sserr, ier = i - 0L, 
     PACKAGE = "biglm")[ 
     c("d", "rbar", "thetab", "sserr")] 

    if(i >= min_window_size){ 
     coef. <- .Fortran(
     "REGCF", np = p, nrbar = nrbar, d = qr$d, rbar = qr$rbar, 
     thetab = qr$thetab, tol = qr$tol, beta = beta., nreq = p, ier = i, 
     PACKAGE = "biglm")[["beta"]] 
     out[i - min_window_size + 1] <- coef. %*% X[i + 1, ] 
    } 
    } 

    out 
} 

# assign function to compare with 
func <- function(formula, data, min_window_size){ 
    sapply(seq(min_window_size, nrow(data)-1), function(x) { 
    mod1 <- lm(formula, data = data[1:x, ]) 
    pred1 <- predict(mod1, newdata = data[x+1, ]) 
    pred1 
    }) 
} 

# show that the two gives the same 
frm <- y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 
r1 <- biglm_wrapper(frm, dat, 25) 
r2 <- func(frm, dat, 25) 
all.equal(r1, r2, check.attributes = FALSE) 
#R> [1] TRUE 

# show run time 
microbenchmark::microbenchmark(
    r1 = biglm_wrapper(frm, dat, 25), r2 = f2(frm, dat, 25), 
    times = 5) 
#R> Unit: milliseconds 
#R> expr   min   lq  mean  median   uq  max neval cld 
#R> r1 9.976505 10.00653 11.85052 10.53157 13.36377 15.37424  5 a 
#R> r2 1095.944410 1098.29661 1122.17101 1098.58264 1113.48833 1204.54306  5 b