Effiziente Methode zum Erstellen rekursiver Schätzungen außerhalb der Stichprobe zur Berechnung von RMSE in R

Question

Ich versuche, die Out-of-Sample-Prognose-Leistung verschiedener OLS-Modelle zu bewerten. Die einfachste Zeitreihenregression sieht wie folgt aus: Y_t = b0 + b1 * Y_t-30 + e_t

Die Armatur Zeitraum für das Modell ist, das 50 sagen lassen, dann lasse ich das Modell laufen mit dem dynlm Paket

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

In meinem aktuellen Code lasse ich einfach den Index laufen bis zum Ende und dann speichere ich den RMSE des entsprechenden Modells. Aber dieser RMSE ist nicht die Out-of-Sample-Vorhersage, und da mein aktueller Code schon ziemlich langsam ist und nicht einmal genau das tut, was ich will, wollte ich Sie fragen, ob Sie einen Vorschlag haben, der Paket, das ich verwenden sollte, um mein Ziel zu erreichen.

Um es zusammenzufassen, möchte ich folgendes tun: vor

1) führen eine rekursive Regression nach einer gewissen Passperiode (erweitert Fenster, NICHT Roll Fenster)

2) erstellen in einem Schritt out-of-Probe prognostiziert

3) berechnen die mittlere quadratische Fehler dieser Prognosen vs. tatsächlichen Beobachtungen das Modell Leistung zu bewerten

ich habe versucht, dies zu tun, so weit mit einem riesigen for-Schleife und der dynlm Paket , aber da ist Ults sind nicht wirklich befriedigend. Jeder Input wird sehr geschätzt, da ich seit einiger Zeit nach Lösungen suche. Ich werde meinen Beispielcode aktualisieren, sobald ich Fortschritte gemacht habe.

# 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

Nun, als derjenige, der die Frage gestellt, ich beitragen möchte, wie ich mein Problem gelöst:

Als ich die einen Schritt voraus Prognosen nur brauchen kann ich das alles andere und das machte wegzuwerfen Code läuft viel schneller. (von 12 Minuten bis ~ 10 Sekunden pro Modell).

Ich erstellte den vollständigen Datenrahmen (einschließlich Lags) selbst und verwendet lm anstelle von dynlm. Der folgende Code gab mir meine gewünschten Ergebnisse (Ich überprüfte die ersten paar Beobachtungen manuell und es scheint zu funktionieren). Der Code wird von hier angepasst: 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) 
       }) 

Für die RMSE Berechnung verwendete ich diese Funktion

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

Danke für die Hinweise auf CrossValidated es sehr geholfen.

Answer 2

Sie können mit den Fortran Funktionen von

Miller, A. J. (1992), um die Rechenzeit ziemlich viel verringern. Algorithmus AS 274: Routinen der kleinsten Quadrate zu ergänzen jene von Gentleman. Journal der Königlichen Statistischen Gesellschaft. Serie C (Angewandte Statistik), 41 (2), 458-478.

können Sie tun, indem Sie mit diesem Code

# 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